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\  Abstract 

A/ 

This  thesis  extends  the  work  of  Leighton  and  Jones 
which  takes  functions  that  satisfy  special  types  of  diff¬ 
erential  equations  and  determines  an  interval  on  which  the 
functions  either  have  zeros  or  attain  bounded  values. 
Theorems  for  locating  zeros  are  proved  for  functions  of  a 
single  variable  and  functions  of  several  variables  with 
illustrative  examples.  The  applications  of  matrix  equations 
to  constrained  optimization  problems  are  described.  An 
algorithm  for  random  search  technique  for  the  general 
optimization  problem  is  presented  with  a  FORTRAN  V  program 
and  test  problems. 


SPECIAL  NONLINEAR 


OPTIMIZATION  TECHNIQUES 

I .  Introduction 

For  tne  past  twenty  years,  considerable  effort  has  been 
expended  on  the  development  of  nonlinear  programming  theory 
and  algorithms  for  solving  nonlinear  programming  problems  . 

Some  of  t^ese  algorithms  have  been  implemented  on  digital 
computers.  It  is  fair  to  say,  however,  that  solving  a  comp¬ 
licated  nonlinear  programming  problem  by  a  computerized  non¬ 
linear  programming  algorithm  is  an  automatic  process.  Unlike 
linear  programming,  where  computerized  algorithms (variations 
on  the  simplex  method,  usually)  have  long  been  able  to  solve 
problems  of  large  size,  nonlinear  programming  is  still  in  its 
infancy  as  regards  its  ability  to  guarantee  solutions  to 
problems  of  even  moderate  size. 

A  major  barrier  for  solving  nonlinear  programming  problems 
is  the  lack  of  a  computationally  oriented  way  of  representing 
nonlinear  functions  of  n  variables .  The  algorithms  are  often 
not  as  efficient  as  they  could  be  because  of  the  inability  to 
compute  automatically  quantities  related  to  the  complicated 
nonlinear  functional  relationships  that  describe  the  models. 

For  example,  the  accurate  and  speedy  computation  of  first 
derivatives  is  a  usual  requirement  for  algorithms  which  solve 
system  of  nonlinear  equations. 
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There  is  a  "basic  dichotomy  in  programming  algorithms: 
they  may  be  designed  to  converge  to  local  or  global  minima. 

A  necessary  condition  for  a  point  to  be  a  local  minimum  of  a 
differentiable  function  subject  to  constraints  is  due  to  Kuhn 
and  Tucker  (1)  and  might  be  considered  the  fundamental  theorem 
of  mathematical  programming. 

Assuming  that  the  problem  to  be  optimized  is  defined  in 
some  way,  the  various  general  methods  of  optimization  can  be 
conveniently  classified  as  follows: 

1.  Analytical  methods:  which  make  use  of  the  cla  •'.cal  tech¬ 
niques  of  differential  calculus  and  the  calculus  c  iriations . 
These  methods  seek  the  extremum  of  a  function  f(X)  by  finding 
the  values  of  X  that  cause  the  derivatives  of  f (X)  with  respect 
to  X  to  vanish.  When  the  extremum  of  f(X)  is  sought  in  the 
presence  of  constraints,  techniques  such  as  Lagrange  multipliers 
and  constrained  variation  are  used.  For  the  application  of 
analytical  methods,  the  problem  to  be  optimized  must  be  desc¬ 
ribed  in  mathematical  terms,  so  that  the  functions  and  variables 
can  be  manipulated  by  known  rules.  For  large,  highly  nonlinear 
problems,  analytical  methods  prove  unsatisfactory. 

2.  Numerical  methods:  which  use  past  information  to  generate 
better  solutions  to  the  optimization  problem  by  means  of  ite¬ 
rative  procedures. 

A  general  summary  of  computer  codes  for  mathematical  prog¬ 
ramming  that  have  been  tested,  documented,  and  are  available 
to  the  public  occurs  in  ^2) .  Somewhat  earlier  there  appeared 
a  collection  of  FORTRAN  listings  of  optimization  codes,  along 
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with  brief  descriptions  of  the  algorithms  and  their  operations 
(3).  The  potential  user  is  left  to  make  his  own  choice  as  to 
which  method  will  best  serve  his  purpose. 

This  thesis  investigates  new  methods  in  nonlinear  opti¬ 
mization  theory.  The  importance  of  this  study  is  to  be  able 
to  use  these  new  methods  in  making  decisions  concerning  modles 
for  wnich  classical  techniques  do  not  provide  sufficient 
information.  The  overall  objective  is  to  address  and  resolve 
the  new  methods  and  apply  them  to  some  special  problems. 

Chapters  2  and  3  deal  with  objective  functions  of  one 
variable  or  several  variables.  Application  of  given  theorems 
provides  a  domain  of  good  starting  points  for  iterative  methods. 

Chapter  4  describes  the  applications  of  matrix  equations 
to  constrained  optimization  problems.  Use  is  made  of  matrix 
equations  to  obtain  solutions  of  certain  classes  of  nonlinear 
equations . 

Chapter  5  gives  an  algorithm  for  random  search  technique 
for  the  optimization  problems  of  continuous  variables  or  integer 
variables . 

Chapter  6  gives  conclusions  and  directions  for  further 
work  in  the  areas  covered  by  the  thesis. 
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II.  One  Dimensional  Unimodal  Objective  Functions 


The  main  objective  of  this  chapter  is  to  obtain  infor¬ 
mation  as  to  the  zeros,  relative  maxima , relative  minima  and 
bounds  for  one  dimensional  objective  functions.  The  objective 
functions  treated  in  this  chapter  are  assumed  differentiable 
and  defined  on  an  interval  [a,b] .  An  objective  function  f(x) 
will  be  required  to  be  a  solution  of  the  differential  equation 
of  the  form: 

[A  (x)  f  (x)  j'+C  (x)f(x)=0  (2.1) 

where 

(i)  A (x J  and  C(x)  are  both  continuous  functions  on  [a,co  ] 
and  A(x)  >  0  on  [a,b] , 

(ii)  f'=df/dx 

Use  will  be  made  throughout  this  chapter  of  a  class  of 
trial  functions  defined  by  definition  2.1. 

Definition  2.1 

A  trial  function  u=u(x,h)  is  a  real  valued  function 
of  x  and  h,  where  h  is  a  parameter,  h>0  and  u  satisfies  these 
conditions : 

(1)  u  >  0 

(2)  u(a,h)  =  u(b,h)  =  0  and 

(3)  u  has  at  least  one  derivative  with  respect  to  x. 

For  example,  on  [0,h]  ,  h>0  ,  one  trial  function  (and 
the  one  used  predominantly  in  this  chapter)  is  of  the  form: 

u  =  x(h-x) 

Other  trial  functions  which  might  be  used  on  [0,h]  are 


u  =  sin  ( n  x/h) 


or 

u  =  xP(h-x)q  ,p,q  si 

Corresponding  trial  functions  may  be  adapted  to  fit 
any  interval  [a,b]  or  h,kh  ,h>0  and  K>1,  to  extend  the 
domain  of  these  trial  functions. 

Associated  with  the  given  function  f(x),  the  trial 
function  u  and  the  differential  equation  (2.1}  will  be  a 
functional  J(x,a,b)  defined  by 

cl  ^ 

J(x,a,b)  =  -  f  (Au 
b 

where , 


-C  u ‘ 


)  dx 


(2.2) 


u  =  du/dx 
x 

Definition  2.2 

The  functional  J(x,a,b)  above  is  a  function  defined 
on  functions. 

For  example, the  functional  J(x,a,b)  is  a  function  defined 
on  both  f (x)  and  the  trial  function  u. 

Lighton  (4)  proved  the  following  theorem  to  show  that, 
whenever  J(x,a,b)  can  be  made  negative  by  varying  the  para¬ 
meter  h,  then  at  least  one  zero  of  f(x)  must  exist  on  the 
interval  [0,h]. 

Theorem  2.1 


If  f (x)  is  a  function  which  satisfies  the  linear  diff¬ 
erential  equation  (2.1),  and  if  u  is  a  trial  function  on 
[0,h]  such  that  J(x,0,h)  <0  for  the  J(x,0,h)  defined  by(2.2) 


then  £ (x)  has  at  least  one  zero  on  [0,h]  . 

The  idea  here  is  to  make  use  of  this  theorem  and  the 
similar  theorems  in  chapters  2,3  in  optimization.  If  f(x) 
is  the  derivative  of  the  objective  function,  then  the  zeros 
of  f(x)  correspond  to  the  extreme  points  of  the  objective 
function.  The  range  [0,h]  for  which  J(x,0,h)  is  negative 
gives  the  range  for  these  extreme  points.  For  the'  unimodal 
functions  it  will  be  useless  to  search  for  the  extreme  points 
outside  the  range  [0,h]  .  In  most  iterative  methods  which 
require  a  starting  point,  the  range  [0,h]  is  the  recommended 
range  of  the  starting  point  instead  of  searching  over  all 
the  domain  of  the  variable  x.  The  choice  of  the  starting 
point  in  [0,h]  will  decrease  the  computational  time  and  the 
number  of  iterations.  The  following  examples  illustrate 
the  application  of  theorem  (2.1)  for  different  functions  that 
satisfy  equation  (2.1). 

Example  2 . 1 

The  function  y  =  sin(x)  is  known  to  be  a  solution  of 
the  differential  equation 

y  » *  +  y  =  0 

A  comparison  with  equation  (2.1)  gives 
A(x)  =  C (x)  =1 
then,  with  the  trial  function 

u  =  x(h-x) 

and  varying  h  until  J(u,0,h)  becomes  negative,  i.e 

h<*  2^ 

J(u,0,h)  =  J  {  (h-2x)  -  (hx-x  )  |  dx  <  0 

0 
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solution  of  this  inequality  gives 


h>VTo 


thus  the  function  y  has  at  least  one  zero  on  (0  A fa). 


NOTE:  the  actual  zero  occurs  at  x=*  =  3.1416 


Example  2.2 


The  function  f  =  P^(x)  is  the  Legendre  polynomial  which 


is  a  solution  of  the  differential  equation 


(l-x2)f"  -  2xf '  +  n(n+l) f  =  ((l-x2)f’)'  +n(n+l)f  =  0 


2 

Now,  >(x)  *  1-x  is  positive  for  all  |x|<  1, 


C(x)  *  n(n+l) 


choose 


u  *  x(h-x)  ,then 


a  a  2  2 

J  -  J|(l  -  x2)  (h  -  2x)2  -  (n+n)(xh  -  x  )  }  dx  <  0 


the  solution  of  this  inequality  gives 


h  >VlO/  (n  +  n  +  4) 


(2.3) 


Table  I  compares  the  values  of  h  calculated  from  equation 


(2.3)  with  the  smallest  positive  zero  of  P  (x)  calculated 

n 


from  equation  (2.4)  of  the  generalized  Rodrigues  formula 


for  diffe^~nt  values  of  n, 


(n)  (n)  2  i 

(d  /dx  )  (x  -  1) 


(2.4) 


n 

2  n! 


’  r* ,  .•'QTi'viTrc 


V**«rj*T**Tin 


TABLE  I 


Comparisons  of  h  with  zeros  of  P  (x) 


.7906 


.6455 


5423 


zeros  of  P  (x) 
n v  ' 


.  7746 


3400 


.  5385 


Example  2.3  Mathews  (5) 

Let  y  be  the  objective  function  given  by  the  solution 
of  the  ordinary  differential  equation  (2.5) 

y '  '  +  xy  =  0  (2.5) 

In  this  example  we  are  to  approximate  the  first  zero  of  Airy's 
function  y(x)  shown  in  Figure  1.  Matching  of  equations  (2.5) 
and(2.1)  gives 

A(x)  =  1  ,  C(x)  =  x 

Let  u  =  x(h-x) 

J(u,0,h)  =  f  |(h-2x)^  -  x(hx  -  x2)  j  dx  <  0 


The  solution  of  this  inequality  gives 

h  >(20)1/3 

1/3 

then  there  is  at  least  one  zero  for  y  on  [0,(20)  ). 

Komkov(6)  developed  theorem  (2.2)  that  is  a  generali¬ 
zation  of  theorem  (2.1).  In  the  general  case  we  will  consider 


the  nonlinear  equation 


[  A(x)V '  ]»  +  C(x)£ (V)  =  0  (2.6) 

where 

(i)  x  C  )  , 

(ii)  A(x)  >0  for  all  xC[a,b]  and  A(x)CC2  [x,  , «»  )  , 

(iii)  C (x)  C  C  [xa,ea)  ,  and 

k 

(iv)  f(V)  is  of  the  form  f(V)  =  v  ,  k  is  a  positive  integer 
and  k  >  1 . 

Theorem  2.2 

2 

Let  u  be  a  C  [x„o<»)  trial  function,  and  G(u)  be  C1(-«* 

+  )such  that  G  |  u  (a)  |  =  Gju(b)|=0,G|u(Xj^>0for 

all  a<x<b.  Let  g(u)  denote  g(u)  =  G'(u),  and  let 
2 

g  (u)/G(u)  be  a  bounded  function  of  x  on  the  interval  [a  ,b ] 
and  m  denotes 

g(u) 

m  =  max  (  -  — )  (2.7) 

xC  [a,b]  4G(u) 

let 

bf  2 

J (u ,a,b)  =  J  (Au1  -  CG  )  dx  <0  (2.8) 

a 

k 

then  any  solution  V(x)  of  equation  (2.6)  ,  with  f(V)  =  V  , 
k  >1  which  satisfies  V(a)  ^  0  will  satisfy  the  inequality 

.  .  m  1/ (k-1) 

|v  (x)  |  <  ( - ) 

k 

on  some  subinterval  of  [  a,b]  . 


The  following  examples  illustrate  the  applications  of 
theorem  (2.2)  for  different  objective  functions. 

Example  2.4  -Komkov^) 

The  Emden’s  equation  (2.9)  occurs  in  astrophysics.  It 
arises  in  the  discussion  of  a  simplified  thermodynamic  model 
of  a  contracting  nebular  cloud. 

y ' *  +(2/x)y*  +yn  =  0  (2.9) 

When  equation  (2.9)  is  reduced  to  the  form  of  equation  (2.6) 
it  gives 

(x2y ' ) '  ♦  x2yn  *  0 

In  the  special  case  which  occurs  in  physics,  n=5,  the  follow¬ 
ing  solution  is  known 

3c  1/2 


where, c  is  an  arbitrary  nonnegative  constant. 

A  complete  solution  to  the  Emden's  equation  is  not 

known. 

Application  of  theorem  (2.2)  allows  us  to  estimate 

how  close  the  solutions  approach  zero  on  a  given  interval. 

Let  the  objective  function  y  be  given  by  equation  (2.10). 

This  function  satisfies  the  differential  equation  (2.6)  with 

A(x)  *  x2  and  C(x)  *  x2 

2 

choose  G(u)  *  4u 


and  u  =  x(h-x) 


The  technique  is  to  use  equation  (2.8)  to  find  the  smallest 
positive  h  for  which  J  is  negative, i.e 


0  >  J 


)  ‘  Au'2 


C  G  )  dx 


-  h/  \  x2  (h-  2x)  2 
0 


4x2 (h-x)  J  dx 


The  solution  of  this  inequality  gives 

h  >  1.8708 

Komkov  (6)  solved  the  problem  for  h  =  2  and  found  that  m  =  4. 
Hence  every  solution  V(x)  of  equation  (2.9)  with  n  =  5  will 
attain  a  value 


|  V (x)  |  <  (4/5) 1/4  =  .945 

on  some  subinterval  of  [ 0 , 2 ] 

Example  2 . 5  -Komkov (6) 

Consider  the  equation 

y”  +  (n  +  1/2  -  x2/4  )  y  =  0  (2.11) 


with  n  =  6. 

A  solution  of  this  equation  is  given  by  the  parabolic  cylinder 

function  y  (x) ,  which  may  be  expressed  in  terms  of  the  Hermite 
o 

polynomial  H^.(x) 


y6(x) 


-(x2/4) 

e 


Vx) 


(2.12) 


where 

H6(x)  =  x6  -  15x4  +  45x2  -  15 

12 


■  /•  ,V'  •  < 

^i*i.  ■  -  tf-  *  _ 


y^(x)  vanishes  at  x  =  .62  and  at  x  =  1.90 
Komkov  (E)  proved  that  every  solution  of  equation  (2.11) 
vanishes  on  [0,1.8  ].  He  used  fixed  h  =  1.8.  If  we  use  h 
as  a  parameter  and  let 


u  =/Tx/h 


2 

choose  G(u)  =  sin  u 


then  g(u)  =  2sin  u.cos  u 
Hence 

2  2  2 
g  (u)/4  =  sinzu.cos  u  £  G(u)  =  sin^u 

then  theorem  (2.2)  can  be  applied.  Comparison  of  equation 

(2.11)  with  equation  (2.6)  gives 

A(x)  =  1, 

C(x)  «  n  ♦  1/2  -  x2/ 4 


J (u ,0,h)  =  J  (  Au 


CG  )  dx 


/|  (  K/h)2  -  (13/2  -  x2/4)sin2(nx/h)Jdx 


solution  of  this  inequality  gives 


7f2/h  -  13h/4  -  h3/16K  <  0 
using  some  experimental  values  of  h  we  have 
J (u ,0 , 1 . 77)  =  0.01941 
J(u,0, 1 . 78)  =  -0.041074 

Thus,  choosing  h  =  1.78  indicates  that  the  function  y  has  at 
least  one  zero  on  [0,1.78]  which  improves  on  Komkov's  result. 
The  improvement  in  the  result  is  due  to  varying  h  as  a  para¬ 
meter  instead  of  choosing  constant  value  for  h  =  1.8. 
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Theorem  2.3 


Let  the  objective  function  V(x),  or  it's  derivative  be 
a  solution  of  the  ordinary  dif ferentictl  equation  (2.13)  on 
xeja.b] 

\  A(x) .D(V) .V'  |  ’  +  2B (x) .D (V) .V '  +  C(x).f(V)  =  0  (2.13) 

where,  (•)  indicates  differentiation  with  respect  to  x. 

If  the  following  conditions  are  true  : 

(1)  the  function  A  is  continuously  differentiable  with  res¬ 
pect  to  x  and  A > 0  on  [a,bl  , 

(2)  the  functions  D  and  f  are  continuously  differentiable 
with  respect  to  V, 

(3)  the  functions  B  and  C  are  continuous  on  [a,b], 

(4)  f  =  0  only  if  V  =  0, 

(5)  Q(x)  is  positive  definite  for  some  continuous  function 
E (x)  on  [a,b]  where 

r A  'B1 

Q(x)  =  (2.14) 

L-B  E  -* 

(6)  J(u,a,b)  <0  with 


then,  either 

(i)  f  and  V  have  zeros  on  [a,b]  ,  or 

(ii)  u  >  D(df/dV)  on  some  subinterval  p,q  -  [a,b]. 
Proof : 

Let  u(x)  be  given  as  a  trial  function  with  u(b)  =  0 

and  u(a)  =0.  If  V  has  a  zero  on  [a, b] then  the  conclusion 
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follows  immediately . 

Suppose  that  V  f  0  on  [a,b]  and  hence  f  f  0  also  on 
[a,b].  Since  u(a)  *  u(b)  =  0 
then 

K  2 

/  A.D.V* . u  (1/f )  dx  =  0  (2.15) 

a 

Also,  since  Q  is  positive  definite,  then,  by  equation  (2.14) 
it  follows  that 


0<  Ju ' -uDV' (1/f)  ,  uj 


A  -B  u' -uDV’ (1/f) 


B  E 


(2.16) 


From  equation  (2.15)  and  equation  (2.16)  the  following  inequ¬ 
ality  holds: 


2  uDV'  2 

uDVir  -  2uB (u '  -  +  Euz 

f  1 

+  (ADV’u2/f)  1  dx 


b/|  Au'2  -  2Au  '  V '  D  (1/f )  +  Au2D2V,2(l/f2)  -  2Buu' 


+  2Bu2DV'(l/f)  +  Eu2  +  (ADV') ’u2(l/f) 


+  ADV 


’(u2/f)'  } 


*  J"  i  Au  '  2  -  2Buu'  +  Eu^  -  2Auu  '  V  '  D  (1/f ) 

a  ' 

+  Au2D2V ' 2 (1/f2)  +  2Bu2DV ' (1/f)  -  2BDV'u2 (1/f) 
+2ADV,uu' (1/f)  -  ADV’u2V’ (df/dV) (1/f2)  -  Cfu2(l/f)|dx 
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dx 


=  /  |  Au’ 2  -  2Buu '  +  (E  -  C)u2  | 

b  2 

+  /  Au2V,2(l/f2) (D  -  D^|-)  dx 


© 


= J  (u ,  a,  b)  +  J  Au2V*  L  ( 1  /£ 2 )  (D2  -  D )  dx 
a  QV 


Since  J (u ,a ,b)  <  0 

tften,  it  follows  that  either 

(i)  V  =  0  and  f  =  0  ,  or 

f\  j  r 

(ii)  D  >D-^-  on  some  subinterval  p,q  <=  [a,b]. 

Example  2.6  -Mahfoud  (7) 

1/3 

Let  V(x)  =  (sin  x)  '  be  a  given  function  that  satisfies 
the  ordinary  differential  equation 


(V2  V') '  +  V3/3  *  0 


(2.17) 


A  comparison  of  equation  (2.17)  with  equation  (2.13)  gives 

A  =  1,  D  =  V2,  C  =  1/3,  f  =  V3  and  B  =  0 

Let  u  =  x(h-x)  be  the  trial  function  for  h>0,  then 

h  2  2 

0  >J(u,0,h)  =  J  (Au'  -  Cu  )  dx 

0 


=  hJ  |  (h  -  2x)2  -  (1/3)  (hx  -  x2)  |  dx 


Solution  of  this  inequality  gives 

1/2 


h  >  (30) 


then,  either 


(i)  V  has  at  least  one  zero  on  [  0,(30)^^2],  or 

(ii)  D2  =  V4>V2  (df/dV)  =  V  .3V2  =  3V4  which  is  impossible. 


16 


i/2 

Then,  V  has  at  least  one  zero  on  [0,(30)  ]. 

Example  2.7 

Let  V(x)  be  an  objective  function  which  satisfies  the 
ordinary  differential  equation 

7  2n+l 

(vV)  •  +  (k/x)V  =  0  (2.18) 

where  n  is  a  positive  integer  and  k  is  a  real  constant. 

A. comparison  of  equation  (2.13)  with  equation  (2.18) 

gives 


2n 

k/x  , 

2n+ 1 

A  =  1  ,  D  = 

V 

,  c  = 

and 

f  =  V 

u  =  x(h  -  x) 

be 

the  trial 

function. 

Then 

0>  J (u,0,h)  = 

h 

5 

0 

(Au'2_ 

2  , 
Cu  ) 

dx 

= 

h 

5 

0 

{  (h  - 

2x)  - 

(k/x) (hx  -  x2)2  |  dx 

Solution  of  this  inequality  gives 

h  >  4/k 

Then,  either 

(i)  V(x)  has  at  least  one  zero  on  [0,4/k  ],  or 

(ii)  D2  =  V4n>  D(df/dV)  =  V4n  (2n  +1)  which  is  impossible 
for  positive  values  of  n. 

Then,  V(x)  has  at  least  one  zero  on  [  0,4/k]  . 

Example  2.8 

Consider  the  Cebysev  polynomial  function  which  is  a 
solution  to  the  Cebysev  ordinary  differential  equation 


where,  n  is  a  positive 

integer . 

Here  we  have 

A  =  1  -  x2,  D  =  1, 

B  =  x/2  and 

C  =  n2 

Let  u  =  x(h  -  x)  be  the 

trial  function 

and  choo 

E  (x)  = 

const-ant 

such  that 

r  A 

-Bq  rl-x2 

-*/2  1 

Q  - 

- 

L-B 

EJ  L -x/2 

E  J 

is  positive  definite. 

Assume  n  =  2  ,  then  to  find  h>  0  such  that  y  has  at  least  one 
zero  on  [0,h],  compute 

h  .  ? 

0  >  J(u, 0,h)  =  /  {  Au,z  -  2B  uu'  +  (E  -  C)u^  |  dx 

0 

h  -  2 

-  /  {  (1  -  xz)(h  -  2x)  -  2(x/2)(h  -  2x)(hx-xZ) 

0  2  2 

+  (E  -  4)  (hx  -  x  )  \  dx 

Solution  of  this  inequality  gives 


h> 


im 


1/2 


172 


(15  -  2E) 


(2.19) 


then,  for  a  suitable  choice  of  E <  15/2  that  makes  Q(x)  posi¬ 
tive  definite,  we  have  either, 

(i)  y(x)  has  at  least  one  zero  on  [0,h  ]  ,  or 

(ii)  D2  =  1  >  D(df/dV)  =  1  which  is  impossible. 


Thus  y  has  at  least  one  zero  on[0,h],  h  is  given  by  (2.19). 


Algorithm  for  Numerical  Solution 

This  algorithm  determines  the  extreme  point  of  a  uni- 
modal  function  (the  zero  of  the  derivative)  numerically. 

The  inputs  of  the  algorithm  are: 

1-  The  derivative  of  the  objective  function. 

2-  The  limits  of  the  variable  x£(a,b)  . 

3-  The  accuracy  required  to  determine  the  extreme  point. 

The  output  of  the  algorithm  is  the  value  of  h  that  make 
the  functional  J  negative.  The  value  of  J  is  calculated  by 
the  Romberg  integration  algorithm. 

Step  1:  Initialization. 

Step  2:  Choose  initial  value  h  =  a  +  d  where  d  is  the  step 
size  for  the  increment  in  h. 

Step  3:  Calculate  J  by  Romberg  algorithm  for  the  limits  a 
and  h . 

Step  4:  If  J  is  negative  go  to  step  8. 

Step  5:  If  J  -  0  put  h  =  h  +  d. 

Step  6:  If  h  -  b  go  to  step  10. 

Step  7:  Go  to  step  3. 

Step  8:  Print  h  and  h-d. 

Step  9:  Go  to  step  11. 

Step  10:  Print  "No  zero  between  a,b. 

Step  11:  Stop. 

The  zero  of  the  derivative  of  the  objective  function 
lies  between  h-d  and  h.  The  accuracy  of  determining  the  zero 
depends  on  d.  The  computer  program  of  the  above  algorithm  is 


written  in  FORTRAN  V  and  is  attached  in  appendix  A. 

The  numerical  solution  of  Example  2.3  by  the  algorithm 
is : 

2.700  ^  h  -2.701 

for  d  =  .001 

The  case  of  multimodal  functions  of  one  variable  is  just 
an  extension  of  the  case  of  unimodal  function.  The  idea  is  the 
transformation  of  the  origin  to  the  first  h  that  makes  J  nega¬ 
tive.  Following  the  same  steps  as  in  the  case  of  unimodal 
functions  but  with  the  new  origin  at  h  the  algorithm  can  locate 
the  second  extreme  point.  Other  extreme  points  can  be  located 
by  further  transformations  of  the  origin. 


III.  Special  Objective  Functions 
of  Several  Variables 

In  this  chapter  objective  function  W(X)  =  W(xl ,x2 , . . ,xn) 
will  be  considered  where  the  variables  xl,x2,..,xn  are  bounded 
The  notion  of  an  n-dimensional  trial  function  u (xl ,x2 , . . ,xn ,hl 
h2,..,hn)  will  be  used  and  is  an  extension  of  the  u(x,h)  used 
in  chapter  2.  The  theorems  are  presented  for  the  general  case 
of  n  variables.  Proofs  will  be  given  only  for  the  case  of  two 
variables  for  simplicity  (see  Figure  2);  however,  they  can  be 
extended  to  .apply  to  functions  of  several  variables. 

3.1  Objective  Functions  that  Satisfy  First  Order  Partial 
Differential  Equations 

The  objective  function  W(X)  is  required  to  be  a  solution 
of  the  first  order  differential  equation 

n 

£  [  Ai(X).B(W).W(X)]xi  +  C(X).f[W(X)]  =  0  (3.1) 

i*i 

where 

X  =  (xl,x2 , . . ,xn) ,  withen  the  region  R  defined  by 
R  (xl  ,x2  , . .  ,xn)  :  O^xl^al ,  0-x2* a2  , . .  . , 0-xn4an  |  ,  (3.2) 

[A]xi  denotes  the  derivative  of  A  with  respect  to  xi 
Theorem  3.1 

Let  W(X)  be  a  given  objective  function  which  satisfies 
the  differential  equation  (3.1)  where 

(1)  Ai(X)  >  0  and  continuously  differentiable  over  R  defined 
by  (3.2),  i  =  1,2, ...,n. 


21 


(2)  B(W)  is  continuously  differentiable  function  of  W. 

(3)  f[W(X)]  =  0  whenever  W(X)  =  0. 

(4)  C(X)  is  continuous  on  R. 

Let 

hi  h2  hn  n  7 

J  =  J  /  ...  J  [  V  Ai.u ^  -  C.G(u)  ]  dX  <  0  (3.3) 

°o  o  fn 

where 

(i)  u  =  u  (xl  ,x2 , . . .' ,xn,hl  ,h2 , . .  .  ,hn)  is  a  trial  function 
differentiable  on  R  and  is  zero  on  the  boundry  of  R,  and 
positive  in  the  interior  of  R. 

(ii)  G(u)  >  0  in  the  interior  of  R  and  vanishes  on  the  bound¬ 
ary  of  R. 

(iii)  dX  =  dxl  dx2  . ...dxn. 
then,  either 

(a)  W(X)  has  at  least  one  zero  on  H,  where 

H  =  |  (xl ,x2 ,  . . .  ,xn)  :  0*xl6hl,  0£x2^h2 , .  . . , O^xn^hn  j  (3.4) 
or 

n 

(b)  Yj  Ai[(l/4)g2(u)B2W2  -  BG(u)Wxi(df/dW  ]  >  0  (3.5) 

i  »•  r 

on  some  subregion  H  -  R,  where 

g (u)  =  dG(u)/du 

Proof: 

To  prove  theorem  (3.1)  we  will  use  a  special  case  for 
n  *  2  and  the  proof  can  be  extended  for  the  case  of  n-variables. 
Assume  W(X)  does  not  vanish  on  R  and  consequently  f[W(X)]  will 
not.  For  simplicity  wesHall  use  A,u,g,B,f,lV  and  G  without 


arguments.  Also  derivatives  will  be  represented  by  subscripts 
(du/dxl  =  ux^) . 

By  Green’s  Divergence  theorem  (8)  we  have 


hlr  h2 


J  J  [ (Al.B.W.G/f )xl  +  (A2.B.W.G/£)x2  ]  dx2  dxl  = 


then,  the  following  inequality  holds: 


h2r  hi/-  2  2 

0  <  J  J  [  Al(uxl  -  gBW/2f)  +  A2(ux2  -  gBW/ 2f ) L 

+  (AIBWG/f)  +  (A2BWG/f )  ]  dxl  dx2 

A  X  a  4 

Ii  2  liX 

=  J  j  [Al.uxl  -  Al.gBWuxl.l/f  +  Al.g2.B2.W2(l/4f2) 

+  A2.u22  -  A2.g.B.W.ux2(l/f)  +  A2.g2.B2.W2(l/4f2) 

+  (Al.B.W) (f .g.uxl  -  G.Wxl.df/dW) (1/f2)  +  (Al.B.W)xl(G/f) 

+  (A2.B.W)(f .g.ux2  -  G.Wx2.df/dW) (1/f2) 

+  (A2.B.W)x2(G/f)*  ]  dxl  dx2 

h2/»  hi..  7  2  A.  7  o?? 

=  J  J  [A1  *uxl  +  A2.ux2  -  C.GJ  +  [  £(Ai/0(g2B  tC/4 


-  B.W.G.VL. .df/dW  )]  dxl  dx2 


J(u,hl,h2)  +  Jf  (Ai/£2)(g2.B2.W2.l/4 


-  B.W.G.W  .df/dW)  ]  dxl  dx2 
Since  J(u,h2,hl)<0  by  hypothesis,  then,  either 

(a)  W(X)  has  at  least  one  zero  on  H,  or 

(b)  22  Ai[g  -B  •  K  /4  -  B.G.WX.. df/dW  ]  >  0 

i»*  t- 

on  some  region  H  -  R 

Theorem  (3.1)  says  that:  if  a  region  H  *=■  R  can  be  found 
for  which  a  trial  function  u  makes  J (u ,hl ,h2)  <  0,  then,  W  will 
either  vanishes  somewhere  on  R  ,  or,  the  inequality  (3.5)  is 
satisfied  on  some  region  of  R. 

3.2  Objective  Functions  that  Satisfy  Elliptic  Partial 


Differential  Equation 


Theorem  (3.1)  applies  to  certain  functions  that  satisfy 
the  first  order  partial  equation  (3.1).  For  other  class  of 
functions  that  satisfy  the  elliptic  partial  differential  equ¬ 
ation  we  have  to  use  the  following  theorem: 

Theorem  3.2 

Let  the  objective  function  V(X)  satisfy  the  elliptic 
partial  differential  equation 


1 1 

^2  [ai;j(X).B(V).Vxi  ]  +  b  (X)  .  f  [  V  (X)  ]  =  0  (3.6) 


wherJ 


(1)  a^j (X)  is  an  element  of  a  symmetric  positive  definite 
matrix  A,  and  has  continuous  partial  derivatives  with  respect 


(2)  X  =  (xl ,x2 ,x3 , . . . ,xn)  is  a  point  belonging  to  a  smooth 
closed  bounded  region  R;  R  is  defined  by  equation  (3.2). 

(3)  B(V)  and  f(V)  are  nonlinear  functions  of  V(X)  such  that 
f(V)  ^  0  when  V  0. 

Let  the  functional  J(u)  associated  with  V(X)  and  equa¬ 
tion  (3.6)  be  of  the  form 


hnr  h.2r  hi/-  7 

J(U)  =  0 }  ■  ■  J  0 /[Sai3CX)UxiUx;  -  b™u  1 

i*3 


dX 


where 


(i)  dX  *  dxl  dx2 . dxn 

(ii)  the  trial  function  u(X)  has  continuous  partial  deriva¬ 
tives  with  respect  to  xl ,x2 , . . . . ,xn  and  it  is  also  zero  on 
the  boundary  of  R  and  positive  in  the  interior  of  R. 

then,  for  J(u) <  0,  either 

(a)  V(X)  has  at  least  one  zero  on  some  region  H  =  R,and  H  is 
defined  by  equation  (3.4),  or 

(b)  B2/f2(V)>  B (df/dV)  for  H  C  R. 

Proof: 


We  will  use  the  same  methodology  as  in  proof  of  theorem 
(2.1)  with  the  same  simplicity  assumptions  for  n  =2  . 

Assume  V(X)  does  not  vanish  on  R,  and  consequently 
f[V(X)]  will  not.  Since  A  is  a  positive  definite  matrix,  let 

all  ai2 
a21  a22_ 


A  = 


C 


then 


1 


Vuxi  »  x2  “*',xi*i/ 

Lux2 


Also  by  Green's  theorem  (8)  we  have 


uBV  .1/f-j 

X1  ] 
uBVx2.l/fJ 


>0 


(3.7) 


h2y*  hi  -  .  _ 

J  J{l“2CanBVxl  .  a12BVx2).l/f]xl 

*  [u2(a21BVxl  ♦  a22BVx2).l/f]x2  }  dxl  dx2  =■  0  (3.8) 

Combination  of  equations  (3.7)  and  (3.8)  gives  the  following 
inequality : 


0<I  -  J  J{* ll“xl  •  2aUuBuxlVxl-l/f  *  *nuWxrl/i2 

*a12Uxlux2  -  a12uBux2Vxr1/f  *  a12UBuxlVx2’1/f 

*  a12u2B2VxlVx2-1/f2  *  a21ux2uxl  -  a21uBux2Vxr1/f 

'  a21aBuxlVx2-1''f  *  a21u2B2vxlVx2-I/f2  *  a22ux2 

-  2a22uBux2Vx2-1/f  *  a22uWx2-1/f2 
.  [u2(anBVxl  -  a12BVx2).l/f 

.  [u2(a21BVxl  .  a22BVx2).l/£  ]  }  dxl  dx2 

h2 *  hi/-  . 

■  J  J  {  (  Z/aijuxiuxj)  *  <-2anuBuxivxr1/f 

27 


.  •  .  .  .  .  k'fc 


■  _>  /»  .*<  A  .V  *.  .%  A  A  aVV  /  ' 


+  ailU  B  Vxi  -  a12uBux2Vxl'1/f  ■  a12uBuxlVx2'1/f 

*  a12u2B2vxlVx2'1/fZ  '  a21uBux2Vxl'1/f  -  a21uBuxlVx2'1/£ 

*  a21u2B2vxlVx2'1/f2  -  2a22uBux2Vx2.l/f 

*  a22u  B  Vx2-1^f  ^ 

*  (aXlBVxl  +  al2BVx2)t2fuUxl  "  “2-df/dV.Vxl)(l/f2) 

+  (a21BVxl  +  a22BVx2)(2£uUx2  •  u2-«/dV.Vx2)(l/f2) 

*  tu2/£)[(auBVxl  *  a12BVx2)xl  ♦  (a^BV^  ♦  .^SV^)  }dxldx2 


and  by  equation  (3.6)  we  have 


h2-  hi 


>  1  =  f  f  [(  V*  a.  u  .u  . 

0 J  0 J  f-f  .  X1  XJ 

\  >Js| 


)  -  bu‘ 


+  V  (a  u  V  V  .1/f  ) 
xj  xi  Xi 
i  ,  j  =1  J 


2  2 

' (B  -  B.df/dV)]  dxl  dx2 


=  J(u,hl,h2) 


h2-  hi/ 


Since 


nir  ,v^  2  7  2 

/  /  a.  Vv  -V  .  (l/fz)  (B  -  B.df/dV)  dxl  dx2 

3 J  oJ  13  X1 

£  a‘- 

A'., 


i,  j  =  l 

.Vx^Vv^  is  positive  because  A  is  positive 


xj 


definite,’ •'and  since  J<0  by  hypothesis,  it  follows  that  either 


(i)  V  has  at  least  one  zero  on  R,  or 

2  2 

(ii)  E  /f  >  B.  df/dV  on  a  subregion  H  ^  R. 

The  following  example  illustrates  how  to  apply  that 
theorem  for  a  practical  function. 

Example  3.1 

Let  V(x,y)  be  a  function  which  satisfies  the  different¬ 


ial  equation: 


v  +  V  v  +  kW  +  =  0 

xx  yy  1  2 


we  have 


b  =  1,  B  =1  and 


f  =  k-jV  +  k2V' 


(3.9) 


choose  the  trial  function: 

u  =  sin  (fix/h^)  sin(ny/h2) 
a  comparison  of  equations  (3.9)  and  (3.6)  gives 

all  =  a22  =  1  and 

a12  -  a2i  *  0 

which  satisfies  that  A  is  a  positive  definite  matrix.  Also, 


h->  h 


\  V  (u2  +  u2  -  U2  )  dx 

J  0J  x  y 


h2_  hl, 


qJ  0  J  {  K*/hl)  cosfrx/hj)  sin(*y/h2)  ] 


+  [  (n/h2)  sin(7ix/h1)  cos(^y/h2)  ] 

-  [sinfxx/hj^)  sinfay/l^)  ] 2  j  dx  dy . 
calculatin  of  the  double  integration  gives 


(h1h2/4)(^2/h2  +  7f2/h2  -  1  ) 


(3.10) 


Since  J^O  and  hp  are  positive,  then 
C7T2/hf )  +  (7T2/h2)  <  1 
By  theorem  (3.2),  either 

(a)  V(x,y)  has  at  least  one  zero  on  H  -  R  defined  by 
H  *  |  (x  ,y) :  0  -  x  -  hp  0  -  y  -  h  |  ,  where ,  hj  and  h? 

satisfy  the  inequality  (3.10),  or 

(b)  B^/f2  >  B,df/dV  which  gives 

l/(k1V  +  k2V3)2>  kx  +  3k2V2 

i .  e . 

V2(k  +  k2V2)2  (kx  +  3k2V2)  <1  C3- 

for  k  =  0,  the  inequality  (3.11)  gives 
3k2V8  <  1 

which  gives  an  upper  bound  of  V 

V  < (3k2  )’1/8 

for  k  =0,  the  inequality  (3.11)  gives 

V2k*  <  1 


which  gives  an  upper  bound 


of  V 

V<(k^ 


) 


-1/2 


If  for  some  values  of  ^  and  k2  the  inequality  (3.11) 
does  not  hold,  this  makes  case  (a)  true  and  there  will  be  at 

least  one  zero  on  H. 


Algorithm  for  Numerical  Solution 

This  algorithm  finds  values  of  hj  and  h2  that  make  the 
functional  J  negative.  Since  J  is  a  double  integral  (for  the 
case  of  two  variables),  that  can  be  defined  by 


Step  1:  Initialization. 

Step  2:  Input:  lower  bounds  A  and  C  and  upper  bounds  and 

U£  for  x  and  y  respectively;  initial  values 
and  h2*,  incremental  values  d^  and  d2  for 
x  and  y  respectively. 

Step  3:  Find  the  approximation  of  J  by  any  multiple  integ¬ 
ration  algorithm  for  the  lower  bounds  A  and  C  and 
the  upper  limits  h.^  and  h2. 

Step  4:  If  J  'y  0,then  let  hi  =  h-^  +  di  and  h2  =  h2  +  d2, 
otherwise,  go  to  Step  6. 

Step  5:  If  hi  4  Ui  and  4  U 2,  go  to  Step  3,  otherwise,  go 
to  Step  7 . 

Step  6:  If  J  <  0,  print  hi  and  h2- 

Step  7:  Stop. 


If  the  value  of  J  will  be  negative  for  x  =  hi  and  y  =  h^, 
then,  the  extreme  point  of  the  function  will  be  in  the  range 


hi  -  di  4  x  4  hi  and 
h2  -  d2  4  y  4  h2 


The  algorithm  is  written  in  FORTRAN  V  and  is  attached  in 
Appendix#  .  The  double  integral  is  approximated  using  the 


Composite  Simpson's  Method. 

For  the  case  of  multimodal  functions,  the  idea  is  just 
the  transformation  of  the  axes  to  the  point  (x,y)  =  (hi,h2) 
where  (hi,h2)  is  the  point  determined  by  the  previous  routine 
The  first  transformation  is  used  to  locate  the  second  extreme 
point,  and  so  on. 

Figure  3  shows  the  first  extreme  point  in  region  A.  The 
origin  is  shifted  to  The  second  extreme  point  is 

located  in  region  B. 


3.3  Use  of  the  Generalized  Inverse  Technique  in  Optimization 
Now  we  consider  a  general  matrix  A  of  order  mxn  and 
rank  lc  which  may  be  less  than  min(m,n)  and  raise  the  question 
whether  an  inverse  exists  in  some  suitable  sense.  This  nat¬ 
urally  depends  on  the  purpose  for  which  such  an  inverse  is 
used. 

If  A  is  a  nonsingular  matrix  of  order  m,  the  solution 
of  the  linear  equation 

AX  «=  B 

where  B  is  a  mxl  column  vector,  is  given  by 

X  =  A_1B 

where  A  *  is  the  inverse  of  A  (i.e.,  AA1  =1).  We  ask  the 
question  whether  a  similar  representation  of  the  solution, 
that  is  ,  of  the  form 

X  =  GB 

is  possible  when  A  is  a  singular  square  or  rectangular  matrix. 

If  there  exists  a  matrix  G  such  that  X  =  GB  is  a  solution  of 
AX  =  B  for  any  B  such  that  AX  =  B  is  a  consistent  equation, 
then  G  does  the  same  job  (or  behaves)  as  the  inverse  of  A, 
hence  may  be  called  a  generalized  inverse  of  A. 

Definition  3 . 1 

Let  A  be  an  mxn  matrix  of  arbitrary  rank.  A  generalized  irw 
of  A  is  an  nxm  matrix  G  such  that  X  =  is  a  solution  of  the 
equation  AX  =  B  for  any  B  which  makes  the  equation  consistent. 
Lemma  3 . 1 
G  exists  iff 


AGA  *  A 


(3.12) 


Rao  and  Mitra  (9)  proved  the  following  theorm: 


Theorem  3.3 

Let  A  be  of  order  mxn  and  G  be  any  generalized  inverse 
of  A.  Further  let  H  =  GA.  Then  the  following  hold; 

(a)  A  general  solution  of  the  homogeneous  equation  AX  =  0 


X  =  (I  -  H)  Z 

where  Z  is  an  arbitrary  vector. 

(b)  A  general  solution  of  a  consistent  nonhomogeneous  equa 
tion 


AX  =  B 


(3.1 


X  =  GB  +  (I  -  H) Z 


(3.1 


(c)  A  necessary  and  sufficient  condition  that  AX  =  B  is 
consistent  is  that 


AGB  «  B 


(3.1 


3.3.1  Common  Zeros  of  Nonlinear  Functions 


The  idea  here  is  to  use  the  generalized  inverse  technique 
in  constrained  optimization  problems.  The  following  two  exam¬ 
ples  illustrate  this  technique. 

Example  3.2 


Given  the  nonlinear  function 
F(x,y)  =  xy2  +  y3  -  x2  =  0 
with  the  constraint 

V  (x  ,y)  =  x  +  y  -  1  =  0 


(3.16) 


(3.17) 


the  problem  is  to  find  common  zeros  of  both  F(x,y)  and  V(x,y) 
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Writing  equation  (3.17)  in  the  matrix  form  given  by 
(2.13)  gives 


where , 


A 


(3.18) 


The  generalized  inverse  of  A  (see (10)  for  calculation  of  G) 
is 


This  can  be  checked  by  applying  lemma  (3.1)  as  follows: 


AGA 

Applying 


(1  ,  1) 


U  ,  1) 


the  consistency  condition  (3.15) 


=  A 
gives 


AGB  =(1,1)  (1)  =  (1)  =  B 

Then  by  theorem  (3.3),  the  general  solution  of  equation (3 . 14) 
is 

X  -  GB  +  (I  -  H) Z 

where 


then,  if  we  assume 


we  have 
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•(:;] 


(l)  + 


5  -  .  5  t  ru 


]Q 


r  .  5  +  .  5u  -  .  5v  -| 

L  .  5  -  .  5u  +  .  5v  J 

A  Taylor  series  expansion  of  F(x,y)  gives 


(3.19) 


F(x,y)  =  F  (0 , 0)  +  [x.F  (0,0)  +  y.F  (0,0)] 

X  / 

+(1/2!) [x2.Fxx(U,0)  +  2xy.Fxy(0,0)  +  Fyy(0,0)  ] 
+(1/3!)[x3.Fxxx(0,0)  +  3x2y.Fxxy(0,0)  +  3xy2 . Fxyy (0, 0) 


+  y  .F  (0,0)  ]  + 

yyy 


=  (1/2) (X  ,  y) 


+  (1/6) (x  ,  y) 


f2 

0  1  rx- 

L  0 

°J  lx 

r° 

6y“|  rx- 

L  o 

6y-lLy. 

(3.20) 


Substituting  in  equation  (3.20)  by  the  general  solution  given 
by  equation  (3.19)  gives 

2 

F(x,y)  =  -(.5  +  .5u  -  .5v)  +  y(.5  -  . 5u  +  .5v)  =  0  (3.21) 

Since  the  general  solution  given  by  equation  (3.14)  is  valid 
for  all  values  of  Z,  [theorem  (3.3) (b)],  then  choosing 


■a  -a 


and  substituting  in  equation  (3.21)  gives 

-  .  25  +  . 5y  =  0  or , 

y  =  .5 

Using  the  constraint  (3.17)  gives 

x  =  .  5 


Then  the  common  zero  of  equations  (3.16)  and  (3.17)  is 


o 


(x,y)  =  (.5, .5) 

The  next  example  illustrates  the  same  technique  but  for 
equations  of  four  variables. 


Example  3.3 


Find  the  common  zeros  of  the  pair  of  nonlinear  functions 
F(x,y,u,v)  =  x2  -  2y2  -  u2  +  w2  =  4  (3.22) 

V(x,y,u,v)  =  x  -  2y  +  u  +w  =8  (3.23) 

Writing  equation  (3.23)  in  the  matrix  form  given  by 
equation  (3.13)  gives 

V  =  (1, -2  ,1,1) (x,y,u,w)T  =  (8) 
where  T  denotes  the  transpose  of  the  matrix. 

Assume  we  know  the  generalized  inverse  G  of  A  to  be 
G  =  (  3/8,  -1/8,  2/8,  1/8  )T  (3.24) 

which  satisfies  lemma  (3.1)  and  the  consistency  condition. 


where 


H  =  GA  = 


Assume 


then 


X  = 


(3.3),  the  general 

solution 

(I  - 

H)  Z 

f  3/8  - 

6/8 

3/8 

3/8  " 

-!/ 

8 

2/8 

-1/8 

-1/8 

2/8  - 

4/8 

2/8 

2/8 

L  1/8  - 

2/8 

1/8 

1/8  . 

z2  , 

23, 

T 

z4) 

3_ 

5/8 

6/8 

-3/8 

-1 

+ 

1/8 

6/8 

1/8 

2 

-2/8 

4/8 

6/8 

1 

-1/8 

2/8 

-1/8 

-3/8 

1/8 

-2/8 

7/8 


zl 

z2 

z3 

z4 


(3.25) 
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.  V. 

-V. 


j 


is 


s> 

V' 


K.- 


fm*  J 
!-*0 


y.v 


i 


P«l 


K7 


►•5 


ft 


& 


F  =  XT 


1 
0 
0 

L  o 


ion  (3.22) 

as 

0 

0 

0  ' 

-2 

0 

0 

0 

-1 

0 

0 

0 

1  . 

for  XT  and 

X 

X  =  (4) 


(3.26) 


(3.26)  gives  the  general  solution  in  terms  of  zl,  z2,  z3,  z4 


Since  the  vector  Z  is  arbitrary,  then  for  a  special  choice  of 

,  T 


Z  =  (0,  0,  0,  0)  ,  the  common  zero  is 
(x,y,u,w)  =  (3, -1,2,1). 


3.3.2  An  Application  in  Iterative  Methods 

One  of  the  best-known  methods  for  solving  a  single 
nonlinear  equation  in  a  single  variable,  say 

f(x)  =  0  (3.27) 

is  Newton's  (also  Newton- Raphson)  method 


Lk+1 


=  x  -  [  f (xR)/f • (xk)  ]  ,k  =  0,1, 


(3.28) 


Under  suitable  conditions  on  the  function  f  and  the  initial 
approximation  xD,  the  sequence  (3.28)  converges  to  a  solution 
of  (3.27);  see,e.g.,  Ortega  and  Rheinboldt  (11)  ,  for  itera¬ 
tive  methods  in  nonlinear  analysis,  and  in  particular,  for 
the  many  variations  and  extensions  of  Newton's  method. 

Newton's  method  for  solving  a  system  of  m  equations  in 
n-variables 

f}(xl,x2 . ,xn)  =  0 

f  2  (xl  ,x2  , . ,xn)  =  0 


(3.29) 


fm(xl,x2 . .  =  0 
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►o 


is  similarly  given,  for  the  case  m  =  n ,  by 

Xk+l  =  xk  '  [f (xk)/f ' (xk)]  ,k  =  0,1,...  (3.30) 

where  f '  (xk)  is  the  derivative  of  f  at  xk,  represented  by  the 
matrix  of  partial  derivatives  (the  Jacobian  matrix) 

f’(xk)  =  dfi(xk)/dxj  (3.31) 

If  the  nonsingularity  of  f'(xk)  cannot  be  assumed  for 
every  x^,  and  in  particular,  if  the  number  of  equations  (3.29) 
is  different  from  the  number  of  unknowns,  then  it  is  natural 
to  inquire  whether  a  generalized  inverse  of  f ' (xk)  can  be 
used  in  (3.30),  still  resulting  in  a  sequence  converging  to 
a  solution  of  (3.29). 

Ben-Israel  (12)  illustrates  the  use  of  generalized 
inverses  in  a  modified  Newton  method  for  solving  the  nonlin¬ 
ear  equations.  Other  applications  of  generalized  inverses 
in  the  iterative  methods  of  nonlinear  analysis  are  in  (12)* 


IV.  Applications  of  Matrix  Equations 
to  Constrained  Optimization 


The  nonlinear  programming  problem  can  be  defined  as 


follows : 


min  f (X) 


subject  to 

g±  CX)  —  0  ,  i  —  1,2,...  ,m<  n, 

hj (X)  ^  0  ,j  =  1,2,. ..,r 
T 

where  X  =  (xl ,x2 ,...., xn)  and  all  functions  f,  g^  and  h^  are 
differentiable.  Note  that  m,  the  number  of  equality  constraints 
(m)  ,  must  be  strictly  less  than  n,  the  number  of  variables. 

If  m  =  n  the  problem  is  termed  overconstrained,  since  there 
are  no  degrees  of  freedom  left  for  optimizing.  Note  that  the 
inequality  sign  in  h^ (X)  -  0  can  be  reversed  by  multiplying 
through  by  -1  without  changing  the  mathematical  statement  of 
the  problem. 

In  this  chapter  we  will  consider  only  the  case  of  equ- 
ality  constraints.  Note  that  the  inequality  constraints  can 
be  changed  to  equality  ones  by  adding  (subtracting)  slack 
(surplus)  variables.  The  solution  first  proposed  by  Lagrange 
in  1760  was  to  form  a  new  unconstrained  problem  by  appending 
the  constraints  to  the  objective  function  with  Lagrange  mul¬ 
tipliers  w^,  i  =  l,2,...,m.  The  new  objective  function 
L(X,w)  is  called  the  Lagrangian.  It  is  defined  on  Em+n, 


which  is  a  higher-dimensional  problem  than  the  original,  since 
it  has  m  +  n  unknowns.  The  price  paid  for  disposing  of  the 


constraints  is  this  higher  problem  dimensionality.  Since 
the  problem  defined  by 


L(X,w)  =  f  (X)  +  £  w  g i  (X) 

i  =  l  1 

is  now  unconstrained,  we  can  apply  the  necessary  conditions 
for  stationarity  given  by: 

m 

dL/dx.  =  df/dx.  +  £  w.  dg./dx.  =0  ,j  =  1,2, ...,n, 

1  1  i=l  1  x  J 

dL/dw^  =  g^(X)  =0  ,i  =  l,2,...,m. 

which  yield  a  set  of  m  +  n  equations  in  m  +  n  unknowns  (X,w) 

*  * 

to  be  solved  for  the  optimal  values  (X  ,w  ) . 

4.1  Riccati- type  Matrix  Equations 

This  section  deals  with  the  solution  of  the  set  of  non¬ 
linear  algebTaic  equations  related  to  the  Riccati-type  matrix 
equations.  The  use  of  matrix  equations  will  avoid  the  itera¬ 
tive  procedures.  To  obtain  these  solutions  we  use  results  of 
Jones  (13)  to  obtain  solutions  of  matrix  equations  of  the 
following  quadratic  form: 

XDX  +  AX  +  XB  +  C  =  0  (4.1) 

in  which  X,A,B,C  and  D  are  nxn  matrices  having  elements 
belonging  to  the  field  C  of  complex  numbers,  and  X  is  unknown, 
Let  R  and  F(R)  be  defined  by 


'-B 

D  " 

ru  m  -| 

,  F (R)  =  ! 

(4.2) 

_-C 

A  . 

LV  N  -J 

where  R  is  2nx2n  matrix  and  U,V,M  and  N  are  polynomials  in 
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A,B,C  and  D.  Let  also  U  be  the  generalized  inverse  of  the 
matrix  U.  The  following  theorems  are  from  Jones  (13). 


Theorem  4.1 

Let  F(q)  be  any  polynomial  of  degree  n  -  1  in  q  with 
coefficients  belonging  to  C  such  that  R  and  F(R)  are  given 
by  (4.2).  Then  a  solution  of 

(X ,  I )  F  (R>)  =  (0,0) 

with  U'l  or  M'l  existing,  or  a  solution  of 

'<] 

with  M  or  N  1  existing  is  also  a  solution  of  (4.1). 

Theorems  (4.2)  and  (4.3)  establish  other  sufficient 
conditions  for  the  existance  of  solutions  of  equation  (4.1 
Theorem  4.2 

Let  R  and  F(R)  be  given  as  in  (4.2)  w'here  M~*  exists 

and 

V  =  VUU  , 

VU  =  NM"1 

then,  X  =  -NM  *  is  a  solution  of  equation  (4.1). 

Theorem  4 . 3 

Let  R  and  F(R)  be  given  as  in  (4.2)  where  M  1  exists 

V  =  NNV  ,  U  =  MNV 

then,  X  =  M  is  a  solution  of  equation  (4.1). 


and 
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Example  4.1 

Solve  the  following  system  of  nonlinear  equations: 


I  - 


-yz  - 

6x  - 

2y  - 

2z  = 

0 

-yw  - 

2y  - 

2w  + 

2x  = 

0 

-wz  + 

2x  - 

2z  - 

2w  = 

0 

-w2  + 

2y  + 

2w  + 

2  z  = 

0 

(4.3) 


This  system  of  equations  can  be  put  in  the  general  form  given 
by  equation  (4.1)  as  follows: 


-3 

-2 


lj 


0 

0 


°1  r 0  0 1 

o  Lo  o  J 


This  gives 


"  3 

-2 

0 

0  " 

R  = 

2 

-1 

0 

-1 

0 

0 

-3 

-2 

0 

0 

2 

1 

Let 


X  = 


[:  i] 


the  characteristic  equation  is 

2 


det (R  -  ql)  =  (q  -  1)  (q  +  1)  =0 

Table  II  shows  that  there  are  3  combinations  of  the 
characteristic  roots  to  be  checked. 

TABLE  II 

Combinations  of  the  chacteristic  roots 

i - 1 - 1 - 1 


root 


(q-D(q-l) 

(q+D(q-l) 


(q+ 1 )  (q  ■  i ) 

(q+l) (q+l) 


Case  1: 


F(q)  =  (q-l)(q-l)  =  q2-2q+l 


F(R)  =  R  - 2R+I 


0 

0 

0 

2 

‘u 

M  ' 

0 

0 

-2 

2 

0 

0 

12 

8 

V 

N 

_0 

0 

-8 

4  . 

. 

_ 

(4.4) 


Since  M'1  exists,  then  by  theorem  (4.1)  the  solution  of 


(X  ,  I)F(R)  =(0,0) 

will  be  a  solution  of  equation  (4.3).  Using  equation  (4.4) 
gives. 


(X  ,  I)F(R)  =  (X  ,  I) 


=(0,0) 


(4.5) 

(4.6) 


U  M' 

L  v  N  J 

which  gives  the  following  pair  of  equations: 

XU  +  V  =  0  , 

XM  +  N  =  0 
Since  equation  (4.4)  gives 

U  =  V  =  0 

then  equation  (4.5)  can  be  satisfied  by  any  solution  of  equa 
tion  (4.6)  given  by 

-10  6 
5  4  J 

Also,  since  N  ±  exists,  then  by  theorem  (4.1),  the  sol¬ 
ution  of 


-  -km-1  .  r 


(4.7) 


,-l 


F(R) 


r «' 

- 1 

o 

_ i 

i — 

i 

X 

- 1 

o 

_ 1 

will  be  a  solution  of  equation  (4.3).  But 
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Ey  following  same  steps  as  in  case  1, gives 

r-io  6] 

x=  L  .4]  (4-10) 

Thus  the  solutions  of  the  system  of  equations  given  by 
equation  (4.3)  are  the  values  of  X  obtained  from  the  three 
cases . 


4.2  The  General  Matrix  Equations 

In  the  case  of  the  Riccati-type  matrix  equations,  it  is 
easy  to  find  all  the  solutions  by  direct  application  of  theorem 
(4.1).  For  the  case  of  the  general  matrix  equations  of  the 
form: 

AX  +  XB  +  C  +  XDX  +  XEXFX  =0  (4.11) 

in  which  X,  A,  B,  C,  D,  E,  and  F  are  nxn  matrices,  it  is  diff¬ 
icult  till  now  to  find  all  the  solutions.  This  section  deals 
with  the  problem  of  finding  only  some  of  the  solutions  of 
the  system  of  equations  given  by  (4.11). 

Theorem  4.2 

Let  F(q)  be  any  polynomial  of  degree  n  -  1  in  q  with 
coefficients  belonging  to  C  such  that  R  and  F(R)  are  given 
by  equation  (4.2).  Then,  X  is  a  solution  of  the  pair  of 


equations : 


DV  -  DNX  -  UEXFX  =  0 


U  =  MX 


(4.12) 


iff  it  is  a  solution  of  equation  (4.11),  provided  that  M 


exists  . 


Proof : 

Let  X  be  a  common  solution  of  equation  (4.12),  F(R)  is 


given  by  equation  (4.2)  and  M“*  exists.  Since  F(R)  is  a 
polynomial  in  R  and  F(R)  commutes  with  R  [see  (14)],  then 


R.F(R)  =  F (R) .R 


that  gives 


■  -  B 

D- 

-u 

M- 

[u 

M  -l 

r-B 

D  ' 

--C 

A  - 

-V 

N  - 

Lv 

N  J 

L-c 

A  - 

(4.13) 

Matching  the  corresponding  elements  in  the  product  of 
the  matrices  gives 

-BU  +  DV  =  -UB  -  MC 
-BM  +  DN  =  UD  +  MA 

(4.14) 

-CU  +  AV  =  -VB  -  NC 
-CM  +  AN  =  VD  +  NA 

Using  the  identities  (4.13)  the  following  is  obtained: 


0  -  DV  -  DNX  -  UEXFX 

=  BU  -  MC  -  UB  -  DNX  -  UEXFX 
=  BMX  -  MC  -  MXB  -  DNX  -  UEXFX 

=  BMX  -  MC  -  MXB  -  (UD  +  MA  +  BM)X  -  UEXFX 

=  -UDX  -  MC  -  MXB  -  MAX  -  UEXFX 
=  -MXDX  -  MC  -  MXB  -  MAX  -  MXEXFX 

*  -M (XDX  +  C  +  XB  +  AX  +  XEXFX) 

Since  M’^  exists,  then 


AX  +  XB  +  C  +  XDX  +  XEXFX  =0  (4.15) 

Starting  with  equation  (4.15)  and  reversing  the  steps 


using 


U  =  MX 
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we  get  the  system  of  equations  given  by  (4.12)  which  com¬ 
pletes  the  proof  of  the  theorem. 


Theorem  4 . 3 


Let  F (q)  be  any  polynomial  of  degree  n^l  in  q  with 
coefficients  belonging  to  C,  such  that  R  and  F (R)  are  given 
by  equation  (4.2)  .  Then  X  is  a  solution  of  the  pair  of  equ8 
tions : 


(4.16) 


-XMC  +  XEXFX  -  NC  =  0  and 
XU  +  V  =  0 

iff  it  is  a  solution  of  equation  (4.11),  provided  that  U'1 
exists . 

Proof : 

Let  X  be  a  solution  of  (4.14),  F (R)  is  given  by  (4.2), 
and  U'1  exists.  Since  F(R)  is  a  polynomial  in  R  and  F(R) 
commutes  with  R,  then  (4.12),  (4.13)  and  (4.14)  hold.  Using 
the  identities  (4.14)  with  (4.16)  gives: 


0  «  - (XM  +  N)C  +  XEXFV 

*  -XUB  -  XMC  +  XUB  -NC  +  XEXFV 
=»  -X  (UB  +  MC)  +  XUB  -  NC  +  XEXFV 
=  -X (BU  -  DV)  +  XUB  -NC  +  XEXFV 
=  -XBU  +  XD(-XU)  +  XUB  -  NC  +  XEXFV 


=  -XDXU  -  AXU 

-  XBU 

+ 

AXU  +  XUB  -  NC 

+  XEXFV 

*  -XDXU  -  AXU 

-  XBU 

- 

AV  -  VB  -  NC  + 

XEXFV 

*  -XDXU  -  AXU 

-  XBU 

- 

CU  +XEXFU 

*  - (XDX  +  AX  + 

XB  + 

C 

+  XEXFX)U 

Since  U  1  exists,  then 
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AX  +  XB  +  C  +  XDX  +  XEXFX  =  0 


(4.17) 


Starting  with  (4.17),  conversing  the  steps  and  using 

XU  +  V  =  0 

we  get  the  system  of  equations  given  by  (4.16)  which  completes 
the  proof  of  the  theorem. 


Example  4.2 


Find  the  common  solutions  of  the  system  of  equations: 


AX  +  XB  +  C  +  XDX  +  XEXFX  =  0 

where 


the  characteristic  function  is 


|  R  -  ql  |  =  (q2  -  3)(q2  -  4) 

=  (q  -  31/2)(q  +  31/2)(q  -  2)  (q  +-2) 

Table  III  gives  the  all  possible  combinations  of  roots 


for  the  polynomials  in  q  of  degree  2. 


-4q  +  4 


TABLE  III 


Let  us  try  one  of  the  combinations.  Choose 


F(q)  =  q  -  4 


then. 


F(R)  =  R  -  41 
0  2 
0  -1 
0  0 
0  0 


0 

3 

0 

2 


3 

0 

0 

-1 


-  [U  M] 

L  U  KT  -I 


,-1 


Since  M"  exists,  applying  theorem  (4.2)  gives 


U  =  MX 

[°  2h°  3] 

L0  -1J  l-3  0J 

which  gives  the  solution 
'0  1/3' 

X  = 

0  1/3 

Since  X  given  by  equation  (4.19)  satisfies 


DV  -  DNX  -  UEXFX  =  0 

then  this  solution  is  a  solution  to  the  given  problem. 

Trying  other  combinations  can  lead  to  some  other 
solutions  of  the  problem.  The  difficulity  here  is  that 
this  technique  does  not  guarantee  to  find  all  solutions 
for  the  general  matrix  equations. 


V.  A  Random  Search  Algorithm 

The  essence  of  this  chapter  is  the  development  of  a  multi¬ 
dimensional  random  search  algorithm.  The  technique  will  solve 
both  integer  and  continuous  nonlinear  optimization  problems. 
There  are  many  optimization  procedures  which  enable  one  to  find 
the  minimum  of  a  unimodal  function  in  n- space.  If  the  function 
is  differentiable,  global  minimum  may  be  obtained  through  the 
use  of  derivatives.  However,  the  problem  of  global  optimization 
of  multimodal  function  has  received  comparatively  little  atten¬ 
tion,  more  so  when  the  function  in  question  is  non-dif ferentia- 
ble . 

As  a  general  principle,  the  accuracy  with  which  a  proced¬ 
ure  locates  optima  improves  with  the  number  of  functional 
evaluations.  In  principle,  however,  one  seeks  a  balance  bet¬ 
ween  a  degree  of  certainty  and  the  cost  of  implementation.  A 
procedure  which  locates  optima  with  great  precision  and  cer¬ 
tainty  would  be  practically  worthless  if  it  requires  economically 
unfeasible  number  of  calculations. 

There  are  several  search  methods  presently  utilized  to 
seek  global  optima.  Borah  (15)  compared  ‘the  method  of  random 
search  with  the  method  of  systematic  search  and  concluded  that 
the  random  search  method  is  better.  His  conclusion  was  based 
on  the  fact  that  the  systematic  search  suffers  from  the  stand 
point  of  being  computationally  uneconomical  in  that  the  number 
of  evaluations  increases  geometrically  with  an  increase  in  the 
number  of  variables.  Also  the  method  of  systematic  search 


does  not,  in  general,  obtain  all  the  local  minima.  This,  in 
turn,  may  lead  to  some  doubt  as  to  where  the  actual  global 
minimum  occurs.  Among  the  random  search  methods,  are  those 
suggested  by  Brooks  (16),  Becker  (17)  and  Price’s  CRS  method 
(18).  The  simple  random  method  accepts  the  optimum  function 
value  as  global  optimum  after  making  a  specified  number  of 
trials  randomly  selected  from  the  domain.  The  stratified 
random  search  method  divides  the  domain  into  a  number  of  sub- 
domains  of  equal  size  and  selects,  at  random,  a  trial  point 
from  each  subdomain  and  each  time  keeps  the  optimal  function 
value.  The  procedure  is  repeated  many  times.  Some  improve¬ 
ment  on  the  simple  random  search  is  provided  by  Becker  (17). 

His  procedure  begins  with  a  simple  random  search  over  the 
domain.  Instead  of  retaining  the  single  point  with  the  opti¬ 
mal  function  value,  he  retains  a  predetermined  number  of  points 
with  optimal  function  values  in  each  trial.  If  the  number  of 
trials  is  sufficiently  high,  the  retained  points  tend  to  cluster 
around  some  optima.  Then  a  mode  seeking  algorithm  is  used  to 
group  the  points  into  discrete  clusters  and  to  define  the 
boundaries  of  the  subregions  each  embracing  a  cluster.  The 
clusters  are  graded,  by  searching  in  each  for  the  retained 
points  with  the  lowest  function  value  and  then  rated  accor¬ 
ding  to  the  relative  values  of  the  cluster  minima.  The  entire 
procedure  is  then  repeated  using  as  the  initial  search  region 
that  subdomain,  defined  by  the  mode  seeking  algorithm  around 
the  best  cluster.  The  user  may  choose  to  examine  also  the 
second  best  cluster,  or  indeed  all  clusters,  according  to  the 


extent  of  his  doubt  as  to  whether  or  not  the  global  minimum 
will  be  found  in  the  subdomain  defined  by  the  best  cluster. 
Price  (18)  suggested  the  controlled  random  search  (CRS)  that 
combines  the  random  search  and  mode-seeking  algorithm  into  a 
single  continuous  process. 

This  chapter  studies  the  problem  of  obtaining  global 
optima  of  the  general  functions  (differentiable  or  not)  of 
several  variables.  The  procedure  begins  with  evaluating  the 
given  function  at  pre-determined  number  of  points  selected 
randomly  over  the  closed  bounded  domain.  Suppose  M  points 
are  selected  randomly  over  the  domain  and  the  function  is 
evaluated  at  each  of  the  M  points.  The  minimal  functional 
value  and  the  point  at  which  the  minimum  occurs  (if  the  prob¬ 
lem  is  one  of  minimization)  are  saved.  This  step  is  carried 
out  N  times.  The  resulting  N  points  will  cluster  around  the 
minima.  An  illustration  of  this  aspect  is  shown  in  Figure  (4). 
Suppose  there  are  many  cluster  points,  then  there  is  a  poss¬ 
ibility  that  around  each  cluster  point,  a  local  minimum  may 
exist.  We  develop  a  single  program  to  find  all  the  cluster 
groups  as  well  as  cluster  points.  Using  a  local  optimization 
routine  gives  the  exact  minimum  of  each  cluster.  Thus  the 
global  minimum  is  obtained  by  simple  comparison  of  the  obtained 
exact  minima. 

If  some  of  minima  lie  near  to  each  other,  this  procedure 
cannot  separate  the  clusters,  because  the  radius  of  the  hyper¬ 
sphere  which  embrace  these  cluster  points  should  be  very  small 
and  therefore  many  points  still  remain  outside  of  any  of  the 


hyperspheres.  These  points  which  are  outside,  give  false 
cluster  groups  and  thereby  increase  the  function  evaluations 
later  tremendously. 

After  separating  the  cluster,  the  next  task  is  to  find 
the  actual  minimum  in  each  cluster  group.  Any  local  optimi¬ 
zing  method  may  be  used.  However,  Nedler  and  Mead  Simplex 
Search  Method  (19)  is  the  most  efficient  one  for  the  nondiff- 
erentiable  functions.  Thi3  simplex  algorithm  requires  n+1 
points  for  n-dimensional  space. 

5.1  The  Nonlinear  Simplex  Method 

The  nonlinear  simplex  method, not  to  be  confused  with  the 
simplex  algorithm  of  linear  programming,  is  a  direct  search, 
descent  method.  It  is  based  on  a  geometric  construct  referred 
to, in  the  case  of  an  n-dimensional  space,  as  an  n-dimensional 
simplex.  The  simplex  method  for  minimization  may  be  summarized 
as  follows: 

1.  Construct  a  regular  simplex  in  the  parameter  space  of  the 
variables  of  optimization  and  evaluate  the  objective  function 
at  each  vertex. 

2.  Find  the  centroid  of  the  simplex  without  the  worst  vertex 
(worst  vertex  being  the  vertex  with  the  highest  objective 
function  value) . 

3.  Define  a  new  simplex  by  eliminating  the  worst  vertex  and 
adding  a  new  vertex  obtained  by  reflecting  the  worst  vertex 
through  the  centroid.  Evaluate  the  objective  function  at  the 
new  vertex  and  return  to  step  2. 

The  elegant  simplicity  of  this  algorithm  has  led  to  a  larg 


number  of  variations  on  the  basic  method.  In  1964  Nelder  and 
Mead  extended  the  simplex  technique  by  allowing  irregular 
simplices  (line  segments  not  necessarily  all  of  the  same  length) 
This  provided  both  a  degree  of  scale  invaring  and  allowed  for 
a  form  of  acceleration  in  the  search.  Their  modifications  to 
the  simplex  rules  include: 

(a)  If  the  reflected  vertex  is  the  best  vertex  (best  vertex 
being  the  vertex  with  the  lowest  objective  function  value) 
then  try  an  expansion  step  to  another  vertex  that  is 
further  along  in  the  same  direction  that  yielded  the 
reflected  vertex. 

(b)  If  the  reflected  vertex  has  the  worst  objective  function 
value,  then  try  another  vertex  that  is  retracted  towards 
the  centroid. 

(c)  If  the  retracted  vertex  has  the  worst  objective  function 
value,  then  contract  the  entire  simplex  by  moving  towards 
the  best  vertex. 

(d)  Stop  the  procedure  when  the  simplex  shrinks  to  a  suffic¬ 
iently  small  size. 


5.2  The  Nonlinear  Integer  Search 

The  development  of  methods  for  integer  variable  progra¬ 
mming  was  initiated  in  the  field  of  linear  programming,  and 
most  work  has  continued  to  the  general  application  within  this 
area.  In  the  field  of  nonlinear  integer  programming  substan¬ 
tially  less  seems  to  have  been  done.  The  majority  of  the  work 
that  have  been  proposed  in  the  literature  to  date  are  centered 


around  one  of  the  following  four  basic  concepts: 

(a)  Rounding-off  the  continuous  optima:  The  most  common 
approach  to  nonlinear  discrete- value  programming  problems 
is  to  treat  the  variables  as  continuous;  then,  once  the 
continuous  optimum  has  been  determined,  an  additional  search 
is  executed  to  find  a  feasible  set  of  discrete-valued  varia¬ 
bles.  The  most  basic  procedure  is  to  select  the  feasible  set 
of  variables  nearest  to  the  optimal  point.  It  is  well  known 
that  this  procedure  can  lead  to  undesirable  answers  and  the 
point  selected  may  not  represent  the  discrete  optimum. 

(b)  Adaptation  of  nonlinear  optimization  techniques:  In 
many  cases,  nonlinear  discrete  valued  programming  problems  are 
considered  as  a  class  of  nonlinear  programming  problems  in 
which  the  discreteness  of  variables  is  one  of  the  restrictions. 
Many  nonlinear  programming  techniques  of  this  type  have  been 
devised,  mostly  for  solving  engineering  design  problems. 

(c)  Linear  approximation  and  binary  representation:  Another 
class  of  approaches  considers  the  discrete  nonlinear  problem 
as  a  primarily  integer  programming  problem  with  nonlinear 
characteristics  that  may  be  linearized.  Some  of  the  approaches 
use  techniques  involving  piecewise  linear  approximation.  Others 
involve  the  transformation  of  a  nonlinear  function  into  a  poly¬ 
nomial  function  of  binary  variables  and  then  the  transformation 
of  the  polynomial  function  into  a  linear  function  of  binary 
variables.  Several  integer  linear  programming  techniques  have 
also  been  devised  for  solving  nonlinear  integer  problems  direc¬ 
tly,  usually  after  conversion  to  binary  variable  problems. 


Some  of  these  techniques,  assume  separable  functions  with 
certain  monotonicity  properties. 

(d)  Direct  search  methods:  Nonlinear  discrete  search  methods 
represent  another  class  of  approaches  which  are  different  from 
all  the  methods  stated  before.  Since  nonlinear  discrete  search 
occurs  only  over  the  set  of  points,  it  requires  the  function 
values  at  these  points  only.  A  reliable  discrete  search  tech¬ 
nique,  however,  is  not  easily  devised  due  to  resolution  ridge 
difficulties  and  no  procedure  demonstrated  to  be  reasonably 
reliable  has  appeared  in  the  literature. 

5 . 3  Problem  Formulation 

The  new  search  technique,  to  be  presented  will  consider 
the  functions  defined  on  a  closed  and  bounded  domain.  The 
constraints  are  of  the  type  which  serve  to  bound  the  individual 
variables.  That  is: 

ai  4  xi  4  bi 

where,  x^  is  the  ith  decision  variable. 

A  unique  feature  of  this  algorithm  is  its  ability  to  handle 
problems  of  either  continuous  variables  or  integer  variables. 
The  algorithm  has,  therefore,  been  named  COIRS,  an  acronym  for 
Continuous  Or  Integer  Random  Search.  Most  of  the  existing 
techniques  for  solving  nonlinear  integer  variables  programming 
problems  were  designed  to  solve  only  some  particular  problems 
with  specific  structures.  The  objective  i  mction  must  be  ex¬ 
pressed  analytically  and,  in  addition,  many  techniques  contain 


several  special  requirements  on  the  nature  of  the  function 
such  as  continuity,  differentiability  or  concavity.  Moreover, 
the  procedures  are  often  exceedingly  complicated  which  make 
them  difficult  to  understand  and  to  program  for  the  computer. 

The  formulation  of  the  optimization  problem  to  be  solved 
by  the  presented  algorithm  requires  the  objective  function  to 
be  expressed  analytically.  The  problem  of  V  variables,  P 
equality  constraints  and  Q  inequality  constraints  could  be 
stated  in  one  of  the  following  forms: 

(a)  For  continuous  variables: 

Minimize:  F(X) 

subject  to:  a^  -  x^  -  b^  ,i  =  1,2,...,V 

G-j(X)  =0  ,j  -  1,2 . P 

(X)  -  0  ,k  =  1,2,  .  .  .  ,Q 

where,  X  =  (x^ ,X2 , . . . ,Xy) 

(b)  For  integer  variables: 

The  same  objective  function  and  constraints  as  in  (a)  are  used 
in  addition  to  the  constraint 

x^  integers,  i  =  1,2,.  .  .  ,V 

If  bounds  on  some  variables  are  not  given  in  the  problem, 
"artificial"  bounds  are  supplied  in  a  way  that  reflects  the 


user's  guess  of  these  bounds. 


5.4  The  COIRS  Algorithm 

Step  1:  Initialization. 

Step  2:  Input  data. 

Step  3:  Generate  predetermined  number  of  random  points  M 
(inside  the  domain)  for  one  random  search. 

Step  4:  Evaluate  the  function  at  all  the  generated  points. 

Step  5:  Find  the  minimum  of  the  evaluated  functions. 

Step  6:  If  the  number  of  the  random  searches  done  is  less 

than  the  predetermined  number  of  random  searches  N, 
go  to  Step  3. 

Step  7:  Find  the  lowest  minimum  of  the  N  minima  determined 
in  step  5. 

Step  8:  Form  a  cluster  from  some  of  the  minima  determined 
in  Step  5  inside  a  hypersphere  with  center  at  the 
lowest  minimum  not  yet  included  in  any  cluster. 

Step  9:  If  the  number  of  points  in  the  clusters  formed  is 
less  than  M,  go  to  Step  8. 

Step  10:  Use  a  local  optimization  technique  to  find  the  exact 
minimum  near  the  center  of  each  cluster. 

Step  11:  Find  the  global  minimum  by  comparing  the  minima  in 
Step  10. 

Step  12:  Stop. 

Remarks : 

(1)  The  generation  of  points  in  Step  5  is  done  by  a  random 

generator  for  both  continuous  and  integer  numbers  inside  the 

domain  of  the  problem. 


(2)  The  local  optimization  technique  used  in  Step  10  is  the 
Nelder  and  Meads  Simplex  technique  that  gives  a  continuous 
minimum.  For  the  case  of  integer  problems  the  technique  is 
modified  to  find  the  lowest  value  of  the  function  at  integer 
hypercube  around  the  continuous  vertices.  The  number  of  verti¬ 
ces  of  the  hypercube  are  2V,  where  V  is  the  number  of  variables 
For  example,  if  the  continuous  vertex  exists  at 

X  =(1.2,  3.7,  4.9) 

the  vertices  of  the  hypercube  will  be 


Xl 

U, 

3, 

4) 

X2 

= 

(2, 

3, 

4) 

X3 

= 

(1, 

3, 

5) 

X4 

= 

(2, 

3, 

5) 

X5 

= 

(1, 

4, 

4) 

X6 

= 

(2, 

4, 

4) 

X7 

= 

(1, 

4, 

5) 

X8 

= 

(2, 

4, 

5) 

The  vertex  with  the  lowest  value  of  the  function  will  enter 
the  simplex  in  the  place  of  the  point  that  has  the  worst  value 
of  the  function. 

(3)  The  constraints  are  supplied  by  the  user  in  a  subroutine 
that  determines  if  the  generated  points  in  Step  3  satisfy  the 
constraints  or  not. 

The  program  is  written  in  FORTRAN  V  and  is  attached  in 


Appendix 


5.5  Testing  the  Algorithm 

The  procedure  of  testing  a  new  search  technique  is  to 
demonstrate  that  the  new  technique  will  solve  the  test  problems 
solved  by  the  established  techniques.  To  demonstrate  the 
ability  of  the  new  technique  one  may  then  present  a  test  prob¬ 
lem  that  can  not  be  solved  by  the  established  techniques  but  can 
be  solved  by  the  new  COIRS  technique,  or  one  may  instead  solve 
the  established  test  problems  using  fewer  objective  function 
evaluations . 

The  following  three  examples  are  from  Borah  (15) : 

Example  5.1 

2  2 

Minimize  f(x,y,z)  =  (x  -  y  +  z)  +  (-x  +  y  +  z) 

+  (x  +  y  -  z)2 

subject  to:  -1  -  x,y,z  -  1 

The  actual  solution  of  this  problem  is: 
f  =  0 

at  x  =  y  =  z  =  0 

The  solution  by  the  COIRS  algorithm  after  1962  functional 
evaluations  is: 

f  =  . 4466x10 " 7 

at  x  =  .000113 

y  =  .0001585 


z 


- .00002563 


Minimize  f(x,y,z)  =  9  -  8x  -6y  -  4z  +  2x2  +  2y2  +  2z2 

+  2xy  +  2xz 

subject  to:  0  6  x,y,z  -  1.5 

The  actual  solution  of  this  problem  is: 

f  =  0 

at  x  =  y  =  z  =  0 

A  comparison  of  the  solution  given  by  COIRS  algorithm 
and  Borah's  solution  (15)  is  shown  in  Table  IV. 


TABLE  IV 

Comparison  of  COIRS  with  Borah's  Results 


COIRS 

Borah 

Number  of  functional  evaluations 

1826 

14144 

f 

. 6537x10  7 

.  953674xlU-4 

X 

1 

1.013671875 

y 

1 

.9921875 

z 

.999 

.986328123 

Example  5.3 

2  2? 

Minimize  f(x,y,z)  =  9  -  8x  -  6y  -  4z  +  2x  +  2y  +  2z 

+  2xy  +  2xz 

subject  to:  0  -  x,y,z  -  1.5, 

x,y,z  integers 


The  actual  integer  solution  of  this  problem  is: 


f  =  0 

at  x  =  v  =  z  =  1 

Using  the  COIRS  algorithm  for  the  case  of  integer  variables, 
the  solution  after  800  functional  evaluations  is: 

f  =  0 

at  x  *  y  =  z  =  | 

Example  5 . 4 

Minimize  f(x,y,z)  =  100 [z- (x+y) 2 / 4 ] ^  +  (1-x)^  +  (1-y)2 
subject  to:  0  -  x,y,z  6  1.5  , 

x,y,z  integers 

The  actual  solution  for  this  problem  is: 
f  =  0 

at  x  =  y  =  z  =  1 

Using  the  COIRS  algorithm  for  the  case  of  integer  variables, 
the  solution  after  540  functional  evaluations  is: 

f  =  0 

at  x  =  y  =  z  «  1 

To  complete  the  testing  of  the  algorithm,  many  other 
problems  have  to  be  solved  and  compared  with  the  solutions 
by  other  algotithms.  This  may  indicate  additional  features 
of  the  COIRS  algorithm. 


VI.  Conclusions  and  Directions 
for  Further  Work 

There  are  several  computer  applications  of  iterative 
techniques  which  will  iterate  towards  zeros  of  a  function. 

The  initial  starting  value  for  these  schemes  must  be  chosen 
carefully.  If  the  initial  point  is  not  close  to  an  actual 
zero,  the  iterative  techniques  often  diverge. 

The  ease  of  performing  the  calculations  given  by  the 
functionals  on  the  trial  functions  makes  these  techniques 
for  determining  whether  a  function  vanishes  or  is  bounded 
in  a  certain  region  attractive.  Other  choices  of  trial 
functions  could  be  made,  perhaps  giving  better  results  in 
certain  situations  than  other  trial  functions. 

The  theorems  and  examples  presented  in  this  work  in 
chapters  2  and  3  can  be  used  to  locate  intervals  on  which 
a  function  has  a  zero.  If  the  given  function  can  be  paired 
with  one  of  the  special  differential  equations  given,  and 
if  a  trial  function  u  is  found  which  makes  the  corresponding 
functional  J  negative  on  a  certain  interval,  then  a  zero  of 
the  function  f  does  exist  on  that  interval.  These  techniques 
have  a  possible  use  in  determining  an  initial  guess  for  a 
starting  point  in  iterative  search  techniques. 

Examples  of  common  functions  from  different  engineering 
sciences  were  chosen  to  demonstrate  the  breadth  of  applications 
possible.  Perhaps  this  method  can  be  applied  to  indicate  the 
existance  of  zeros  of  other  functions  before  valuable  computer 
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tine  is  wasted  to  find  these  zeros.  In  fact,  the  techniques 
obtained  can  thenselves  be  implemented  using  numerical  methods 
on  a  computer  of  reasonable  size. 

There  is  no  doubt  about  the  results  of  techniques  presen¬ 
ted  in  chapters  2  and  3  for  the  unimodal  functions.  For  the 
case  of  multimodal  functions  further  studies  and  analysis  is 
required.  Application  of  given  theorems  for  the  case  of 
multimodal  functions  will  result  in  dividing  the  region  of  the 
function  into  subregions  with  at  least  one  zero  on  each  sub- 
region.  A  comparison  of  this  technique  with  the  known  search 
techniques  is  an  area  for  further  study. 

Analysis  presented  in  chapters  2  and  3  is  limited  to  some 
functions  that  satisfy  special  classes  of  homogeneous  differ¬ 
ential  equations.  The  case  of  nonhomogeneous  differential 
equations  and  other  classes  of  differential  equations  needs 
further  analysis  and  effort.  Other  remaining  aspect  in  this 
area  is  to  compare  our  techniques  with  other  approaches  that 
can  be  used  to  solve  the  problem  of  finding  good  initial  point 
for  iterative  methods. 

For  the  case  of  constrained  optimization,  the  technique 
of  generalized  inverse  was  given  with  illustrative  examples 
in  chapter  3.  The  generalized  inverse  deals  with  the  special 
case  when  the  problem  is  ill  conditioned.  This  technique  is 
applied  to  one  of  the  best-known  methods  for  solving  a  system 
of  nonlinear  equations  in  n-variables.  The  problem  appears 
when  the  matrix  of  partial  derivatives  (the  Jacobian  matrix) 
is  singular. 
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Applications  of  matrix  equations  to  constrained  opti¬ 
mization  problems  is  presented  in  chapter  4.  For  the  special 
case  of  the  Riccati-type  matrix  equations  there  is  no  problem 
for  finding  all  the  solutions  of  the  problem  by  applying  the 
technique  of  matrix  equations.  The  difficulity  appears  in 
the  case  of  the  general  matrix  equations  where  only  some  of 
the  solutions  can  be  obtained  by  this  technique.  The  problem 
of  finding  all  the  solutions  in  the  general  case  is  a  wide 
area  for  further  study  and  analysis.  Some  theorems  have  to 
be  developed  to  guarantee  finding  all  the  solutions. 

The  general  optimization  problem  for  continuous  or  integer 
variables  is  covered  in  chapter  5.  The  COIRS  algorithm  is 
tested  by  some  examples.  Still,  it  needs  further  testing  with 
different  problems.  The  case  of  mixed  integer  problems  using 
the  COIRS  algorithm  is  a  recommended  area  of  further  study  and 
research.  An  algorithm  that  solves  the  general  mixed  integer 
problem  without  any  restrictions  on  the  objective  function  will 
be  a  superior  one  in  optimization. 


Appendix  A 

Computer  Print  Out  to  Find  ji 


PF  ."!(  *S  RllMNRG 

» > ************ x.*x  ****** ».*•**  *#x*.*x»  •,  »•*'** *  x *****  ****** *X 
T.-Tb  PF'OF Rf* F  ]n!  S  T!=r  Vft,.UE  OF  H  *1  W1-.  C.~  THE  F i  sFTIP  if. 

J  JP  NEGATIVE .  The  Pf CULfT  3  ON  P"  :  .  Pr.pr  PV  RQr  •■  PR; 

KL I 

f  =  IDWfcf:  LIP  IT  OF  V 
i:  =  upper  x 

h  =  F'APft*  FTEF 

j  -  tncrpk:£nt  f:p  -■ 

F  X)  =  THE  LNTEGPf.NI  OF  J 

x:  *•**:*•  v  **  *:4  **:*  *:*.*:*  <  !  3  ''  l  ».*:>;  *». x: ».*(:*.'*’!  *»>  ffi.s  **.**;; *  *  ^  *  * 

DINER  ill  ‘fs  T  ( i  2  *  10 
KMftX--6 


6  =  100  * 

H=.l 

D=*001 

0=2 

1 0 L C *■  ’  (H) 

SL'H-  •:  r  ( h  •  f ;  4  f  (  H ,  K  >  1  n . 

HMl=h-5 

IF  (0.11)40, 1C, 8 
DO  V  1=1 ,rthl 

S'JM-aUN+F  ( H .  A+F  LCAT  (()*!..; 

T(l,l)=5Uh*C 
K:  20  K=  2 ,  KMA'X 
OC/2. 

*=M  r 
OUh-O. 

DO  II  1-:  ,M,2 

HUM=Sur.+F  ( H ,  fi  +  F  LOftT  <1)4-0 

T  ( K , )  )  =  T  N - 1 , 1 )  i‘2 , 1  SUM*:C 
RJLIRJ=1  . 

DO  12  J=2*N 
F0UF-.J=PC’JRJ*4. 

GAVE  DIFFERENCES  FOR  LATE'.  1:0.0. C-  FA’  IDF 
T  ( i\  - 1 , ,  I  -  J  )=T(N. 1  (t;-l  ..J-  1  ’ 

’  '  t  .  J)-T  C  ,.1-iHT  (:  -1  .  j-1 iFODPJ  1 .  '• 
j  r  (  T  (  K  ,  J  I  .  L  !  ,  0 . 00002v>5  »  ■  if’  T0  1 
CPF  T  IN! IE 

l  •  J  **'  ♦  X  ‘  2. 

]  F  (  K  ,  (it  •  *.■ )  bC  TO  4  ■ 

IE  H  .  ’  T  ,  V  i  H-  H4 1* 

GO  TC  :! 

FKTNT*.  • 


Appendix  B 

Computer  Print  Out  to  Find  and  H 


PROHRM  HULIN r 

F.’E^L  X  i  ’  i  A  r  t  *  t!  1  i  L‘»  »  ^  *  Hi.’  r  C »  i 
BOMfiON  HI  r * 

UL-- 10. 

U2- 1 0 . 

H  ’  = .  1 
H2-  .  1 
Hi  -  ,  O'. 

D2-.0C' 

A= 0, 

1  (:=C » 

n=  r> 

M-w 

NN=2*N+1 

h*-2*H-1 

hM=0. 
ftF.  =0. 

ao->:  . 

LO  10  1-1  ,N’N 

X  =  n+  <1  —  1  >  j*  H 

hx=(h?-(.  )/r>*K  > 

BN--F  (Xf  C)+t-  (Xji-’2  ■ 

PE=0. 

110=0, 

DO  20  J- ] T  HH 

y=c+j*iix 

Z  =  F  <  X  t  Y 

IK  J.EQ.2K  J/2) >  THEN 
HE=I<E+Z 
ELSE 
B0=B04/. 

EH  Li  IF 
»  CONTINUE 

til  =  i  BN+2*B-  +4*1-  *)#H>'/3 
IE  ( I , EG .  1 .  OF.  •  I . Er< .  i'N )  THE? 
A W-ftN  ,  5 

el  1 1: 

r  (i . EH. 2*1 1  /2)  >  THE '<■ 
ao=ao+a:. 

el  a 

K-r.i:  -fi1 0! 
r  i,r  i< 

END  IF 

10  LCniNiiL 


Appendix  C 


Computer  Print  Out  for  COIRS  Algorithm 


FRCORA-i  'JO IKS 

*.**>-  * >x>  *  *x>  x »  :i*  x  :ts.^  *»x.»».*.x  >  <  st  »  >»».x :  *  « *  x  :<  x  * xx*  x  ».*>  x  x  *  n 
THIS  f-T'T'yvA-;  1;^  'riE  H.MImjm  ,  a  RjHC  T I  -  if-  USIf-'C  f.  T  :*• 
SEA*-.-:-  ■  E.  -  Nit  .  FJNET.-N  '■*:  EE  T  J  P  F  r  EE  r-'Tl  AF  _L 

NOT  , 


x.x  x  x »i$ix*¥S  %*%$.},  x:*»»***~*» * *jm  * x  * x».x 
M-'CNEI'  n  A, :  3 » 2v  » I'J  ■ » pH  AX  •’ IT'.  -1H'  -  ''E;. 
\2C  ■  V,*KAX(20,20>  »>'  (20 > ,  -  2;  : .  r? 
C t c. A >  ( 2 ’  •  10 *3)  »F A  A >'■  1  ( 20 ?  > r'  ■  ■  • \  f  > 
£ L  1  ( 2 v ,  :•  0 )  •  F Z L  ( ? 0  I  .  IF r L  <  2 0 . 3  :  .  C "  A X 


xi.t.*.:*  xx*:**:,  -  J* » x  X  x»>  x 

* 5)  f  F  '  7:0  r 2f’ > . 

<"i>-  '•  •  '■■■ » 20 )  t 

T  2  U  ;  *  2  F  A:  X  '  1  , ? 

2  ::o,3/,Fr,Ax:  J2C  , 


71  2-1  r  3  > » D  ( 6 ) , E  ( 3 )  •  A 1 ' 3  -  2  0 . 7  ' 

I  NTL7  E.T  ;  E  ( 2  >  f '  A  ( 3 , 20  ,  2  0 ) ,  T I:  <6  r  C  •  AX  \  2'j  1 2 


iLX.X'XX  t.xxx’t  *x  XXXxcX.XXx  »:xix  *;m.x  :  <  ***:  ».)!,>  s * a 

ThE  DA'  A  APE  IKE  O'PER  AK'l  LCVEF  T.j'JMiE  0r  THE  VAR] 

#:**  *.  **».*#:  )f:  a.  ?•.##  ».  X*  *:**»*.#>  *#x. #  *  *  *  *.  >' :  *  *  X  X  x  ‘ :  y  *:  a:  *  *  *  ;t, *  x  •  » *  *.  * .  X  * ; 


DAT  A  D/C .  0 » V .  7 , 0 . 0 , 1 .  b » ; .0,1.3/ 

OF  fcN  ;  3 1  FT_L-  '  lFLAlil '  5 
RcW  N 1 1  3 

********* M  ***%**tXMt***  xxx*  **»».<  *.%!  t  *x*x  %%** 

TO  F  •  ND  THE  INTEGER  LIH  '  E  OF  THE  VAP2AKEE 
1 2*- THE  NUMBER  OF  VARIABLES 


11=12X2: 

*1f*tt*xx*.*  #.**  *■'*  x**x:*'*!t  **>  X 


D.  .•  i'_  K  •**  i  *  o 

I F  ( I . . G E  ». i  1 J ~ 3  >  i  GU  TO  1  2: 

IF(  AMD ( IN  U  r  i , 0 )  .CiT  .  0 . 0 ;  12  ■:  K)  =  INT  i  DU  )  .  -  : 
:ii<k)~:n"7  (hk)  ; 

:  i=6 

f  Ai  jl1  =  • :! 

KOuHT -0 

#***#.*#* )M* XX  *»  X» X x 
j  FLA3=1  FOR  THE  INTE3EF:  VAR]  A:-uES 

=  2  FOP  THE  r->'T IK'JOUR  VARIABLES 
M~TH..  NOME r R  OF  KANfOM  OcARCH' 

N=THt  n..  "Rr_R  L'~  f  I  NT  c.  ]  <•*  EC'c  RY  EEAr  -' 
*•»».* tn%% *:» *.**i*»x«*.»x*» x*.* i>  i^niu > >• 

1 F  L  m  1>=  2. 

M  =  lv 
N=1  (• 

X  X  X  *.  X  X  >: »  X  X  *Xx  X  *  jS  *  *  t..»  hi  t 

i>r  r. rF  -  T.  i.  T/F.'  ; N^T  EE. 

»  *  4  /  -*  *>  xxx  * xxx x  ■  **  i  j  » .  1 


■  30  I-  , N *  1 
i  13  k=i  , 

CALL  UNF3F;A<L  ,XR.K  J 
IF  ( IF  LAG . EG .  J  )  GL.  Tu  700 
C-:  TO  y 

:>o  a ( r  ?  i » j > - 1 :>r 

IF  ♦L'  *  I  !•  <  f  ■ '  IA<l-...I,J'--]A(K,I,2>-»: 

A  i ' t .  *  J  ,  *_■ )  * i  A  ( t f  j.  ,  v> ) "  '  *  0 
GO  TO  13 
7  au.  ,j,  i 
3  GOKTIKi;:u 

>r#iM'**********!U'#M.#***)U#M.**  **•**•**#  hu 
CALCULATE  T:=€  FUNC  IQHAL  VALuE' 

KOUNT=::ntJNTER  FOR  FUNCTIONAL  EVALUATIONS 
)i  **##**  *  *  **t* 

L  Nj-J 

K  1  = } 

IF (1  FLAG. £0. 1  )  F ( I f  J  ) -  CAKAI-Kl  ,N1 l 
:  r  (  J  L  *. 0 » 1 0  *  C  F  ( i  ,  J  3  =CAl  <  A ,  Ki  ,  Nl  ' 

KOii^ !  -KOUNT  +  1 

IF (IF LAG. FO.I )  GO  1 0  6 

UP I T F  ( 3 ,  7  A )  ( A  ( K ,  I , J ) ,  !\  =  1 , 12 ) ,  (  (  J  ,  J  » 

26  F ORAAT  ( 4  v ,  3  ( E 1 0 . 4 , 4x  r  4X » E 3  0 . 4/  > 

GO  TO  30 

A  WRIT  E ( 3 , 27 ) ( ] A ( K  ,  I ,  ,1 ) ,  T  =  1 ,121  ,F< J,J> 

27  FQRMT < 4 X .3(317, 3X>,3>, E '  0 . 4 / 

70  28  K=1 , 12 

23  A  ( K ,  I ,  J  >  =  I A  <  t ,  j  ,  J )  *1  *  C 
30  CONTINUE 

**x*  *»'*** »  x%x/*:t*/k*i**x*%*  *%**% ***->.». *x  *#*##*#*.»  **;*%*.x:$i  *.t  *  * 
CALL  SUBROUTINE  TO  FINE  hIN  OF  T*E  10  p 01  NTS  GUST  ECAl.1. 1 M E 
*  *  **.'**  *  X  *  t.  *  *  *  #:  *  *  *  %  %  *  *  %  *  %  V  t:  *:  )j:  %  *  *:  *  t  %  t  *  *:  %  *  *  *  *  *:  *  *  *  * T  *#*'#;#  > * »  v 

K=  J 
T3-?  2 

CALL  F  « I N 1  ( A , t , FnA - . L  A  A  X , K , N , 1 3 ) 

* 4*  >**:#**'* #  **: X *» ***:*'  J*;*  * X' *'#'**: %  **  >**'*:*** * 

CALCULA  :  ir-v-  to  S' ArT  NFVT  RAMKi-  GEARS G 

*#.#*#  *?»  fr**:****  M  #####.»*#*'»  ***•»  ***** * 

a  j=j 


A.;- A.;/ 1C. 

DO  30  t\~* .13,1 

f r  (  IFl  AG ,  LG .  1  )  1 A  i  K ,  1  , 1  l  =  1  A  i  C .  ]  .  J  -  AJ»  ( ••  l)*'*  J 

if <  i Flag. E  > .j)  a <k,i  . i -Ad:. i ,  j)  +nj».  :-i >x*.j 
7t  CONTINUE 


1540= 

1  06 

CONTINUE 

1550= 

GO  TO  200 

1560= 

-  1  1 

K2=K 

1570= 

d  u  i  ii“.  t  ~ . ,  c  j.  , .! 

1580= 

5 LI  ( J ,  I  ■  =  ’-  CL  ( I ) 

1590= 

FCL(i)=0. 

1*00= 

DO  115  111=1,17,1 

1610= 

cl  (  j,:  ,  li  >  =cr  cm  ,ki  ) 

1  620= 

:  15 

C-  CL  ( I , i; 1  =  0 . 

3  670= 

1 1 1 

COM! 1 NUt 

1640= 

or.  TU  2  ’0 

1650= 

121 

K3=K 

1660  = 

DO  126  1=1,63,1 

1670= 

CL1  (.1,1 )  =  =  ()_  ( I  > 

1680= 

FCL ' 1 2  =0, 

-  :Qm- 

DO  125  K .i  =  l,  12,3 

1  ?;'0= 

CU Jr  -  K 1  1  =(:FCL  >7,511 

17-  ;•  = 

-  cr 
.•  /  U 

CD Cl ( 1 , K 1 1=0, 

1720  = 

1  26 

CONTI  Mu; 

1730  = 

Gf t  TO  200 

1740= 

3.31 

M=c 

3  750= 

DO  136  1=1, K 4 , 3 

1760= 

CU  (7,1  > =K CL (  I  > 

1770= 

FClU  >  =C . 

1750= 

HD  1 35  53=1,12,1 

1790= 

CL(J,I,KI)=CFCL( l , K 1 ) 

1800  = 

1  35 

d  CU:  ,  K 1  1=0. 

1810  = 

136 

CONTINUE 

1820= 

200 

CONTINUE 

1630=5 

184Q=C 

UC.7TF.  THE  HINCTION’-OAI  ML  AND  TH‘~  COOF.HNM 

1850=5 

CLUSTER 

1R60=C 

1870= 

220 

I F ( K 1 , EG . 0 )  GO  TO  275 

1880= 

WRI 1 E ( 3 , 222  HU 

1890= 

922 

FORMAT <1 OX U CLUSTER  OF  ',131 

1900= 

HU  224  .1  =  i  ,1.1,1 

1910= 

224 

WRITE  ( 3 ,  *  1  CL1  ( 1 , 1 ) ,  ( CL(  1  ♦  1 »!'.]  ),KI  =  1 .12) 

1920= 

225 

IF (K2.EG.0)  CU  TU  230 

1930= 

WR1  Tr  (3,226) f. 2 

1940= 

226 

FORMAT (10S  'CLUSTER.  UF  '  ,13) 

1950  = 

DO  228  1  =  1 »\2 » 1 

199.0- 

2  9 

WR:TE(L,4)  CU  (2,1  •, (CL I  , K I), K  1  =  1,1 2  . 

1 970  = 

230 

IHK3.cQ»0>  GO  TO  235 

1  990= 

UR I  ' 6(3, 232 If. 6 

1970= 

'j  t  '> 

F riR'16  T  ( :  OX ,  '  Cl  U3  !  EE  (.IF  '  ,  J 

U  00= 

V  234  ! = i ,K3, 1 

701  0= 

234 

Whj:i.-:<3,»  ■  (3.)  (3.  ,  ,l\.  ,U  =  1  .  • 

2020= 

235 

IF  <  K4  ,  Et1 . 0  1  GO  T):  240 

77 


:*C3  <  > 

Wi  :iTE  \  3  t  2' 3  c  ' 

i-  A 

1_!  *.*  A  ? u  • 

236 

F  ORMA'.  •..  m, ' 

CLUSTER  Of",  13,' 

Il/w**'" 

DO  23b’  j=l,K 

4,3 

L*Oc'- 

27H 

WRITE (3, *  'Ll 

1  K  £  r  !  *  (  Lr  L  •  A  Y  •  f 

)'7r) :“' 

240 

Kk-( 

20hO= 

DC  350  2=1,4,:, 

r  1 A  r? r'  •*• 

>33  TC  (302,306,310.314  •  ,i 

2100= 

O  V/ 

iraa.Eo.c;  go  to  350 

22  J  i  0- 

J 1  =  K 1 

21 :  v-- 

I=.J 

-i  n/.- 

I Count =o 

2140- 

4.07 

CALL  SMl-'i.EX  ( LL1  ,  CL ,  •'  YMAX ,  YrftX ,  b ,  I ,  Jl ,  1. 2 , 70;  INI ,  I 

2150=C 

#»#♦*♦***.»  *tt*tt*  ************  x  ***%)**:%  *  >.  i. » 

2160=0 

I COUNT =CQUNTER  FOR  THE  SlHRLE*  ITERATIONS 

?1?0=C 

a  ******  ** 

2100= 

WRITE (3.305) IC0UM 

2.190= 

305 

FORMAT  ( 2>  •  'PK0C''  ITERATION' ,2)  ,  I-*  1 

2700= 

KK  =  FK  +  IL  OiJNI 

2210= 

GO  TC  350 

2220= 

306 

IF  (K2.Eo.C-  >  GO  TO  350 

2230= 

J  1=1(7 

2240= 

I  =  J 

2220  = 

GO  TO  303 

2260= 

3i0 

IF (K3.EO.O )  GC  TO  350 

a.  \7 

J  1=K3 

2280= 

I=J 

2290= 

GO  TO  303 

2300= 

3  j  4 

IE (K4.E0.0)  GO  TO  350 

2310= 

Jl=f(4 

2320= 

1=J 

2330= 

GO  TO  303 

7340= 

350 

CONTINUE 

2350= 

NNT=KK+KOUNT 

2360= 

WRITE(3,355)KNT 

2370= 

355 

FORMAT  (lOX^TOTA!  ITERATION  ,2X,I5) 

2380= 

ST  OR 

2390= 

END 

2400= 

SUXROtJT  1  NE  MIN2 ( CMAX3 , F MAX? , FAX 2 , CF AX2 ,  L . L  ,  1 1 ) 

241 0=C 

#**:*##************** # **#* t, ***: **** ** *******  ** 

2420=C 

THIS  StlRRC'JTINE  FINS  THE  MINIMUM  OF  THE  M  MINIMA 

7430=7 

*****»:**.'  *:**  it*  »****>■  ****  *  ***  ***:*:  ##*.*#**  ».*****♦  * 

7440  = 

DIMENSION  F  A X 2 ( 2 0  > » C m A X ?  v 2 0 , 3 ) ,i MAX3 ( 20 ) ,CF  Ax2 ' 2 

245('= 

11  =  3 

2460= 

DO  10  1  =  1  .K,l 

7470  = 

IF  (FMAXUU  '  .Gl  .FMA*3<  1  )  )  GO  T’’  ]( 

2^80= 

TErR=fhA;-:'3i  -  ) 

249''= 

F  a  A -  ?(!)=■  1 1AT4  ( I ) 

2700= 

FMAX3<1  )  =  7EM! 

251 0= 

0  0  5  f .  =  l  ,  4  *  1 

2520  = 


2  “560= 
2570= C 
256D=C 
2590= 
2.600= 
2610= 


2660  = 
2670= 
268 '*=0 
2690-  C 
2?00=C 
2710= 
2720= 
2730= 
2740= 
2750= 
2760= 
2770= 
2780  = 
2790= 
2800= 
2810= 
2820= 
2830=C 
2*840=0 
2850=0 
2860  = 
2870= 
2880= 
2890= 
2900= 
29 1G=C 
2920=0 


296'.'= 
29-  0  - 
2980= 


Tti1(K)=CW<a»X3(l,K) 

CMAX3(1,K>=CMAX3<1,K1 
CMAX3  !  I  .K  '■  =TEM(  K  ■ 

5  CQ^TIMUF 
10  CONTINUE 

******  Jk'JMiX#*#  *******  **  x* 

STOKE.  THE  CURF.tN"  i  OWEK  VALUE 

*  *  X  X  X  X  X  *  X  a  *  X x  X  »  *  ii.i*  » 4.  **:*#XX* 

^ AX2 l L >=-MA> 3( 1  ) 

PH  AX3  ( 1 )  =  9 . 1 1-  -r "  0 
DO  12  K=i ,3,1 
C  F  AX2( L ,N ) =CMA>3  < 1 , 1  * 

12  CMAX3 ( 1 • K ) =9 , 1E+1 C 
RETURN 
END 

SUBROUTINE  FMIK  ■  :i.,u!,FMA.Xl  »C:M.-V;  ir<.,Ki,u3) 

*  itc  X  if  X 1  It  X  x .  X  *  >v  #  *  .< -  x  X  X  X  x  »:  s»: *  X  X  *  *  X  *  x  *  *.  X  X  * )!:  *  *  *  *  *  x*s:X 
THIS  SUBROU'INt  FINDS  THE  MIN] HUH,  !-F  N  POINTS 
#*x*xx##xx***xx*#x***x**:**x***#*xx*'/-****#**ir>-: 

DIMENSION  K 3 , 20 , 20 )»>-'(  20 . 20 >» FMAX i  (  20  •• , CM  AX 2  <:  2C  ,?  ,TEf, 
L3=3 

i'D  10  j  =•  i  ,  N 1 , 1 

I '•  ( H  ■  j  » L )  ♦  GE  ♦  K 1  1 ,  l  1  5  30  Tc  1 0 
TtM-=H<l,L> 

H ( 1 ,  L ) =H ( I , L  3 
H(  I ,  _)=TEMF' 

DO  5  Kr  =  3 ,1.3,1 
TEH (KK )=D ( KK , 1 » L ) 

&<NK,3 , L ) =D( Kl , , I , L ) 

5  t ( KK , I , L  =TEM ( KF ) 

10  CONTINUE 

»XX»  *X***»!X*»]Mr***fXX)|  XX*  X  XXXX  XX  X*XXX  *  X*  tMM 

STURi  THE  CURRENT  LOWER  VALUE  AMD  COORDINATES 

XX*XXX»XXXXXXXXXXXXXXXXX*X**X*XXXXXXXXX*XXXX* 

FriAXl  <  L  >=H  (!,!.> 

Til  12  KK  =  i,L3,l 
12  CHAX1(L,KK)  =  MKK,1,L) 

RE i URN 
END 

If  tint**  »  ****#**»*XX  *.****  ***X»X**  »  *  X  XXfXXXTXXX  **X  ,U  *  *M  ■< 

LOCA!  O’-TIH.^ATION  .  HIND  Ni'  LD  AMD  M  ADS  S^'L.A.  MET'-.  D 

*****.#**•»  xn  *»#»**)l.****xx**x*x.»:*****s-x»***»*»x*»>xx  v***;* 

SL'BROl'T I f.'r  SHELL  X ;  ALL1 ,  EACl  3  *?‘ZHAX  tZHAX,DI'ti.,A»  4  ■ 

DIMENSION  ACL  1(20.1  OI.CAllI  IS , 20, 3 : ,C2MA A ( 21 , 20 , 3 ; , 
XZMA;  (20, 2C  ;,TI  M~'  ,  IS-;  »  .DD'E*  :  3  ) 

*  ,  X2L  4, 4  >  ,surui)  .PIT  ‘  ,  21  S'  ■  ,C>  ;4,  ■  > .  SR  • . 

X  1  J  MO  1  4  *  A1  ■  ,  f  l.R  (  A  .  ,  t~  ~M  *1  J  •  i‘  (.i'J  A  .  »  r  -  *  (  '  »  v.  L  1  4  ,  . *  ,  *  X  1  L  =  '  ,  £  v 


JO 10=  1  FLAG-2 

3020=  RADD=,9 

3030=  I 4=3 

3040=  ITEF;=0. 

3030=  N2=I4+1 

3060=  HO  3  1=1 rN2» 1 

3070=  BST  ( L. ,  I )  =+ ,  9  1 2=  1 0 

3080=  BO  3  K  =  i  ,  1 4 , 1 

3090=  3  CRGT  ( L » 1 , K 1  =  (  t i  )  !(<> ■; *0 . 9 1  r  *- 1  ■ 

3100=1.  ***** * *********. *****  *********m****» '  *»#.**<«».***»** 

3110=C  FIND  THL  LOWEST  -•CUs-T  AND  IF’-  COORD  nV'FS  BST  (I  .1- 
3  j  20=C  ##******x*********x-***'************#**»*»  ***•#***##* 

3130=  DO  10  1  =  1, 

3140=  IF ( ACL1 (L  , I )  ,Ety»0»91E  +  10)  NO  TO  20 

3*  50=  IF( BST  <  L» 1 ) »Lt * ACLl ( L • I ) )G0  T(,  i ;» 

3160=  T£MP=BST(L,1) 

3170=  BST<L,1)=AClHL,I) 

3180=  ACL1 (L» I )=TEMP 

3190=  DO  8  i’>  =  1 , 1 4 , 1 

3200=  TEMPO  ( K )  =CBST  (  l  ,  1 ,  K  > 

3210=  CKST ; L , 1 , K ) =C ACL 1 < L , T , K ) 

3220=  B  CACLl<t,I,K)=TEMPO(K> 

3230=  10  CONTINUE 

3240=  WRITE(3»*'PST<L,1) 

3250=0  *m**M*:#*#t**#*it:#* 

32&0=C  FIND  THE  SECOND  BEST 

3270=C  Xt*****##*#*:*###*  ******** 

3280=  DO  20  1=1, M,1 

3290=  IF ( ACL1 (L, I ) ♦ ED , F :  1 : T ( L ,  1 ) )  GO  TO  20 

3300=  IF  ( AC-L1  (L, I ) « E Q . 0 , 91E+10 )  BO  TO  20 

3310=  IF (BST(L,2) .  LE*ACL1 (L, I ) )  BO  M  20 

3320=  TEMP=BST<L,2) 

3330=  BST  <L,2 )=ACL1 (L , I ) 

3340=  ACL1 (L , I )  =  T  EMF 

3350=  DO  16  K= 1,1 4,1 

3360=  TEMPO (K)=CPS7 (L,2,K) 

3370=  CBST (L,2,i,)=CA0Ll  (L,1,K) 

3300=  16  CACL1(L,1,K)=TLMP0(K) 

3390=  20  CONTINUE 

3400=  WRITE (3,*) BST (L,2) 

341 0=C  X******  ****  <* *4: ***:** X****** 

3420=C  FIND  THE  THIRD  LOWEST  F0IN": 

343o=C  #**#*XXXX** ».**»***»***  ***** 

3440=  DO  2e.  7  =  1  ,rM 

3450=  IF ( ACL1 L ,  I  - . F  0 ,  i-’RT  t  ,  1  < )  00  TO  2t- 

3460=  IF vACLl <L,I ) .LO.B  -7 < 1 ,2- )  bS  ID  26 

34  ->0=  I F  ( ACL  1 « L ,  i  ).•  it ,  0 . 9 1 Z*  1  0 )  00  TO  „  . 

3430=.  ir(tST(L, 3),:.T  .A1..HL  ,1  '  >  00  TO 

"*490=  TRHpr-DB7  O.  *3; 


(J  cJ 


ACL1  (Lil  E!i'- 
1-C  2  b  K^x  ,34,1 
Th'«;j0(K)*CK-T(L,3,t>  1 
C  B  o  T  ( L.  ,  ,  K ) =  C  A  x  L 1  _  » I  f  f'  } 

CACL1  <i_»  J  >K)=TEMr  0(K  1 
CONTINUE 

URi Th  ( 3  T  i: )  VST  ( L ■  3  .' 

'*'**%*%  %H%%%1t.  %  *  it*.****)  it:  *.  *•  x  t  %  s  **:*»*'  V  *  X  #  *  »:  * 
FIND  THE  FOURTH  LUuH  M  POINT  AT  THIS  1 1  ML 
%  %  **.**:  %  H.%  X;  %  X  *  ***»  x  *:*  *  x  *  *S:  *  *  *xfc  *  i  X  x  ***:  X  **■  >1: 
DO  2?  1=1 ,M,1 


IFTACLl (L, 

1) 

,E0, 

,  BST 

( L , ; ) 

) 

on 

",  0 

±  '> 

IF <  ACL1 (L, 

I) 

.Ely, 

,  DOT 

<  L  ,  2 ) 

) 

uc 

TO 

29 

IF  ( ACL1  <  l.  , 

I) 

»EU . 

•  FT 

<  L  ,  3 ) 

> 

CD 

7C 

OO 

IF ( ACL1 (L , 

1) 

•  DE. 

,  DOT 

( !_  *  4  ) 

) 

CO 

TD 

2  * 

IF  (ALL  I  i.  L, 

I) 

,EU , 

» 0 , 9 1.  E  +  i  0 

) 

BO 

TO 

2^ 

TEMR=BST(L,4) 

BST  <  =  ADL1  ( L»  .•  '> 

AClI<L,I)=TE^' 

DO  28  K= 1,3 4,1 
TEMPO (N  ‘=CH0T(L,4,!\: 

CDGT  (  L  ,  4 ,  p  )  =CiciOul  T  i. » 1  »K  > 

CAClI  1 L ,  i , N )  =  !  EMPO i  ) 

CONTINUE 

WKITi.  (3,*>B?T(u,4) 

*4* #*»•* ****#*iU-*fX**iM jhM* *-*4*»'*».4»-**: »:***»  4 x****** ****>'  4 
THE  SlMf-'LX  FORMED  BY  BST  ( L ,  1 ) ,  BST  ( L ,  2 )  ,  BST  ( i. ,  3 )  t  BET  ’  L  -  *  ■ 
WHICH  I!>  THE  BIGBES:  ONE  SHOULD  BE  REMODEL 
#**#x**X*#****fX**#****#*X*XXx*#X#X»-xx**».)!:y#X##XXM:##!U* 
DO  3b  r .  —  1 » 1 4 » T 
SUM1'  N)=0. 

DC  30  J=J ,N2* 1 

SUM  i  In  )  =SDM  v  K )  +(-£  3T  <  L  ,  I ,  K  > 

DO  36  K= 1,1 4,1 
Xl(K)=SUM(K)/4. 

Dir=o. 

DO  40  I=1.N2,1 

DO  3Y  02-1,14,1 

Y 2  ( I ,  M' )  =C  DOT  <  L ,  I ,  T. 2  •  - X 1  '  K 2  ) 

DIP=Dir+X2U  ,».2)**-2 
CONTINu!. 

Dlc  D-  S  :jF:T  <■  DIF  ) 

I  T  E  K  •  :  '  Sri 

,F  (I.JPfc.LE  .RADD  GO  TU  500 

W1UT  '  3,43 1  B I F 2 

FORMAT  (4a.,  'DIFFERENCE-  '  ,812,5) 

****».»■■  **.#*.****» *  '  %k*i  sHO) r  **,<  *  **x.  *.*.;*:**>. 

vEMOVt  IMF  y  'fiHESY  PR.tM  THE  S]M<\.E>: 

T  HE  H p.HP ,  f  n  1 5  THE  : ,  K  l  ,  *  ) 


y  f ; :::  T 

4000=  L 
4010  = 
4-.  20  = 
4230= 
404  j  = 

4  020  - 
■406  ••- 

"/  1  .  ‘  V*  ~  •_ 

4080- C- 
40?0=C 
4  1  00= 
4’  4  0= 
43.2‘"'=C 
413:'i=C 
4540=2 
4  5  50= 
4'.  60= 
4170  = 
4'.&0=-C 
41 vo=: 
4200  =  C 

/t  '  '  { I 

v  *.  iV 

4220= 


1  Hi:  MEMA'*  ;.u  i  '■'JiO.H  •  *  5  ;  f -H'- 7  ■  l  7  .  F'.  ■  T . 

XXXX>.XX*X»XXX*  »»*  *»:i  X  X  *  xx»:x*  *XX*.\  *  x,  .•  it  *  :*.*»:«  *  *  x  x»4 
PD  43  K  =  1,14,1 
SL'r  2 '  !, )  =0  • 

N’;5=N'"-i 

Tij  4.4  ,  =1 .  N3S  .1 

SU.-12U  )-3:;jH2(!w-C;Lc  ru,?  ,•'■■) 

C^'  1  L  *  K  )  =  SuHt2  ■  P  >/«.-♦ 

xx*  it.  at  xxxxxxx  t  ?.  4  &  4  *  *  i » *  s,* » tx*  *  ■  2 x :  -  ■  • 

F-NIi  7  HE.  IMAGE  OF  P;,'-  ;  L  .  4  0  IHDL'HK  , 

**»:#: 2:* **x.**X*:*0  x *4  H  n* ***  *  *»•>  ,-t  *x*>; * 
i  r;  20  k  =■  1 » :!  4 .  .4 

p  1HC-  ‘  L  ,P  >  =  L  <  HZ?  ■: !. ,  K )  -C  H'-  T  ■:  l  ,4,.  , 

XX*»  XXXXXX  *2  XXX*  XXX*  .XXX  XXX*  *  XX.  XX  />  A  X  <V.  A  *.;!.**  X.t 
C HE C  I!  EXCEEDS  T •-  1  POfMIw  2  -  T-.-  FL :  :  T  IHr. 
xxxx  *xx*  ,  xxxx  j  -xxx  xx  xx  **'*  XX*  **:»  *  *  .s  xu  x  a  x  xx* 

PD  32  f- 1,1 4.1 

I r  (plHGvL , K)  ,L1  ,M  O  '  ,i.5=HP-l- ' 

I F  <  P 1 M •;  L » K ;  ♦  GT » 5  2  -2  +3  >  1  f  Z  4G  -  .  i-.  '  =  -  '  <  F 1 3  / 

>  X X *  X X X  i >.  XX»  X X X X *  *  x: x X X xxxx  >,  ». x X  X  *  ;  X X  v  •< 

EVALUATE  THE  F  L'i.'CT  ]  OH  AT  THP&i  -•■  IP'S 

*  XXXX*  *  XXX  XXX*  XXXX  i  X  X'XXXXXXX  XXX*  X  *  XX* 


4230= 

Ff-T1S(L)=XFC7  ‘  PJh 

G,M  ,L3) 

4240= 

KOL’T  *  K0JT4 1 

,/t ' )  c;  - 

WPITf  13,54  tFPHH-1 

1  ) 

4260= 

54 

FORMAT  MX,  '  FF‘Mfj= 

"  ?  E  3  2 , 4 

4270= 

It  iPPMCIL) , LT . PS 

■  C  1  •  \  nr-  T  r» 

•  ■  M  »  .-  Uw.'  1  U 

4280  = 

]F  <F>-MB<L).L  7  .BE 

•7  (..2)  5  G! 

.1  7  0 

4220= 

IF  (KPMC-'U.PT.BS 

T ( L , 4  '  )  cr 

:  T0 

4300= 

U  &  v  w  K  1‘  it  1  ^  1 1 

/’,3tr‘= 

56 

(Zi  !  (  L  ,  i\  )  =  \  3*VP  (  L  • 

K > -CPST  <  L , 

,4.3, 

7  *  0  -  •  ■ 

XX  X  X  X  X  X  X  XXX  X  s  X*  x 

*  xXXX*.x 

4330=8 

CHECK  I'"  EX'  Et  It- 

DOHA  I  . 

4340=3 

X  *  X  X  X  ■»  X  X  x  X  X  X  X  *  X  x  X  X  X  X  X  X  X 

4  350  = 

DC  58  K=3 - K4 , 1 

4360= 

1 r  ( CF:  ( L ,  K  3  ,27,  PP 

(1  -  CPU, 

,  »  ')  “ 

4370  = 

5G 

II  (CPf  L.K  '  .NT.I.ri 

(F>3)  •  CP' 

:  L  -.- 1- 

43P0=C 

XXPiX.S  xxx*  xx*xvx*xxxxy 

43?  '=C 

tw,  j-r;-  -he  f..ik 

rzTIf.:N 

*4.;>0s.(z 


*  *  1  *  t.f *  t  it  *  X  XXX*  XXX  X*'  • 


*1:1  '*~L 

^520-C 

*  r.  •  ■;  *-•.  r 
t  <.j  w*  v~  v 

«?-.4.j=- 


07  I..  4‘,)(i 
nr  .v'j  14,: 

c* ■  ( c i  ’  \  l  t  k  ■  ■*  i  r-:  ; 

>  * ».»  £*>'**  **:x  »:s  a  i M  Ji 
c:  c  T  r  :  .t  l  r  i  '  ~  -  », 

■Ji.L  il  i  _■  a  L.  !  I'1 

i>H(  ju ***)M  **-*.**  v.  » i  X* 

r:.'  7c  n-t,i4.„ 


4550= 

IF (CwiL  ,!n. LI  no  nn  ■  ci 

4060= 

70 

ir  <c:,'l,io . gt . do (f.+; *. ' 

TT-- 

.« 1=4 

4550  = 

* 

409.:  = 

n.'LM  •  =  ;m < C  r L3 ) 

4600  = 

KG':  1=1. "27  +  3 

A 

o 

vi 

b0  C(L»4  .i=FL'Wiu 

4620= 

or  75  r= i.i 4.1 

4630  = 

“7  cr 

CIS!  fL.  4.  K>=C4 > 

>).440= 

35  TC  400 

4650= 

•i  00 

B?7  ( L .  4  )  =+  F'KGIL  > 

4660  = 

DO  110  1  =1,3.4,  • 

4670= 

lie 

CfeSKl  ,4.K)=F'THG(!..,K  ' 

'4+50= 

GO  TO  4';  ' 

20  v 

Dr  220  K--3  ,  J  4, 1 

220 

tX(L,  :=7.0*C  (L,M-2. 

47; 0=0 

1  >  n:*):*  #*#*  *  **  t  *  S 

4720=C 

Sit  IF  KXCUEDS  DOMAIN 

4730= C 

V*'*»  *X  *  **#•*:#.* 

4740  = 

10  250  K5=l ,14,1 

4750= 

If  ((■>( l. , Kb / » l  7  « DD (K5)  > 

4760= 

250 

I!  (fcX'  L,l\5)  .G  f .  DO  i  f  K .  •  3 

4770= 

+•1  =  74 

4730= 

:.3=! 

4790= 

FED*' < L.  )=XF  CT (EX,  +  ‘  .1.7.) 

4  000  = 

K01IT=K0:)T+1 

40  r  o= 

IK  <■  Ft  X  ( L ) T.l.n(L,3)) 

4820= 

i3i;)TL!I.2;l>  FT*  ID 

4  930= 

.  r 

F*'NM4T(  ;0>.  •  'FF  X='  .  F  •  0 . 

4840  = 

GO  Ti  100 

4850=  .33  0  F’-.T  <L  >-FfX<L  > 

406  0=  DO  200  !."•  .1*34,1 

4k7(>=  T'O  Of  =  ) 

4.?; 0=  40f  :n{ir  .,Af  .E'i.  1 )  GO  7r:  33 
T'</=-  C"  700 

■'0 ’j  —  20:'  DO  3 "•  ■'  1=1.3 


bVr  :L ,  lfO=C 
B:  1  ,  j  , :  )  =  T 


KSTC  f3  ':=V:£T(L  ,4) 

BET  (1  ,4 

i':J  4  l>  -  1  »  4  ,  .1 

CBST <Lj3»\)='.-BST  (i.<rh ) 

C’  4, K)  =  ^Enc.jr:? 

•  r  (itfe.lt,  io.  •)  !••:.  to  v, 

iTts-r* 

PET  <  L  ,  4  >  =  (  F":  T  (  L  ,  D-iBb';  a.  BMP 
GO  TP  400 

ZPi'F  ( L  r  •  >-  l-3T  ( L » '  ) 

DO  50?  '  =1,14,1 

CZMAX  a  ,  1 , K ) =i;BST  P.,3,  K  > 

Uil;l  TP  '3,510)  2H<^>  <!.-3  •• ,  ■' C/ s‘  * 
FOK^TU  0>.  fE  10,4. 7(0.1  0.4, 4>  ) , 
F:s  Ti.  i  N 
END 

F  i ’NT.  i  J  On1  Xi  '  T  ( f. ,  .i. ,  >Ii-  ) 
MOFN-.TOfw  f (4,3) 

1 1-3 

/  r  r  t  -  4  , 0  - '  ,C(CPrV  )~.0»:  •  f 
**;»4 r» , o*t: ■’  ;f  , :> >#* j-l  r : i- , .?  ■  ».* 
p.oxl  I?', :  ■*■:.(  ii-  ,7) 

p:-  t  .j-.-. 

FrJIi 
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